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The  Shearing  Contribution  in  Transformation 
I  Toughening  of  Brittle  Materials 


P.  Sofronis  and  R.  M.  McMeeking 
Department  of  Theoretical  and  Applied  Mechanics 
University  of  Illinois  at  Urbana-Champaign 

March  1985 


SUMMARY 


A  simple  constitutive  law  is  proposed  for  the  description  of  a  ceramic 
composite  which  undergoes  stress  induced  martensitic  transformation.  This  law 
is  used  in  finite  element  calculations  to  investigate  the  shear  effect  on  the 
transformation  zone  near  a  crack  tip.  A  formula  describing  the  stress 
Intensity  factor  change  due  to  the  shear  contribution  of  the  transformation  is 
given.  Significant  loss  of  toughness  is  observed  in  the  case  of  a  stationary 
crack  and  is  attributed  entirely  to  the  shear  component  of  the  transfor¬ 
mation.  On  the  contrary,  the  dilatant  part  brings  about  no  change.  As  the 
crack  grows,  the  wake  of  the  transformed  material  left  behind  the  crack 
constitutes  a  source  of  toughening.  This  toughening  is  due  to  both  dilatancy 
and  shear  in  the  phase  change  and  rises  to  a  maximum  level  just  after  a 
propagation  comparable  with  the  zone  height.  Finally,  it  is  shown  that  the 
shear  component  can  be  Important  when  prediction  of  the  fracture  toughness  of 
the  transformation  toughened  ceramics  are  made. 


1 .  INTRODUCTION 


Transformation  toughening  Is  one  of  the  mechanisms  available  to  overcome 
the  Inherent  brittleness  of  ceramics.  It  is  a  phenomenon  applicable  to 
ceramic  matrices  In  which  Zirconla  (Zr02>  and  perhaps  some  other  materials  can 
be  Incorporated.  To  date  it  has  been  studied  quite  thoroughly  {1-13). 

An  optimally-fabricated  [13]  partially  stabilized  zirconla  (PSZ)  Is  a 
two-phase  ceramic.  Its  microstructure  [5]  at  room  temperature  consists  of:  a^ 
cubic  matrix  which  is  a  high  solute  content  Zirconla  'alloy'  containing  one  of 
the  stabilizers  MgO  CaO,  Y2O3  or  any  of  the  rare  earth  oxides:  and  fine 
coherent  metastable  tetragonal  precipitates  of  low  solute  content  Zirconla 
phase  inside  the  cubic  matrix.  The  stress  induced  martensitic  transformation 
of  those  metastably  retained  tetragonal  particles  to  monoclinic  symmetry  in 
the  stress  field  of  a  crack  tip  is  the  mechanism  responsible  for  the  enhanced 
toughness  observed  experimentally  in  PSZ.  The  phenomenon  of  transformation 
toughening  is  also  observed  in  the  AI2O3  and  ZrC>2  system.  The  pure  tetrogonal 
Zr02  particles  are  retained  metastably  in  the  AI2O3  matrix  by  pressure. 

There  are  two  methods  of  analysis  with  regard  to  the  phenomenon  of 
toughening.  The  first  incorporates  the  energy  changes  accompanying  the 
transformation.  The  second  concerns  the  stress  intensity  factor  (SIF)  changes 
that  take  place  during  transformation. 

Quantitative  analysis  from  a  continuum  mechanics  viewpoint,  based  only  on 
the  dilatational  component  of  the  transformation  with  regard  to  the  SIF 
reduction  was  done  first  by  McMeeking  and  Evans  [19]  and  Budiansky,  Hutchinson 
and  Lambropoulos  [20].  The  analysis  of  the  latter  workers  is  based  on  a 
constitutive  relation  between  the  mean  stress  and  the  dilatation  for  the 
composite  ceramic.  In  both  works  mentioned,  the  predicted  toughening  is 
comparable  with  the  experimental  data,  however,  it  underestimates  them. 
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Recently,  Lambropoulos  [21]  has  suggested  a  promising  constitutive  law 
for  the  composite  including  both  parts  of  the  transformation,  dllatant  and 
shear.  The  treatment  of  the  shear  effect  takes  into  consideration  particle 
size  and  orientation.  The  results  predicted  for  spherical  particles  are  quite 
well  in  agreement  with  experimental  observations  but  are  based  on  a 
transformation  zone  size  and  shape  estimated  from  the  standard  crack  tip 
singular  elastic  field.  That  is,  he  assumes  that  the  zone  shape  is  the  same 
as  regions  in  the  unperturbed  elastic  solution  in  which  the  transformation 
criterion  is  met  or  exceeded.  Changes  in  stress  due  to  the  transformation  are 
not  addressed  in  detail. 

In  this  paper,  we  shall  study  the  shear  effect's  influence  on  the 
enhanced  toughness  through  numerical  calculations  by  means  of  a  simple  stress- 
strain  relation  for  the  composite.  The  model  is  based  on  a  constitutive 
relation  along  the  lines  introduced  by  Budiansky  et  al.  [20]  and  the  condition 
for  transformation  is  dominated  by  the  dilatational  component.  When  the 
hydrostatic  stress  reaches  a  critical  level,  the  transformation  takes  place. 
This  model  is  not  entirely  satisfactory  as  the  shear  effect  is  bound  to 
influence  the  critical  state  for  transformation.  However,  we  regard  this 
paper  with  the  critical  state  determined  solely  by  the  hydrostatic  stress  as  a 
first  step  towards  a  complete  theory.  The  analysis  is  carried  out  around  the 
crack  tip  of  a  long  crack  which  is  imbedded  in  a  composite  material  rich  in 
tetragonal  particles  whose  presence  in  the  matrix  is  defined  by  a  volume 
concentration  Vf.  As  such,  the  composite  can  be  modelled  as  a  continuum  of 
transformable  material. 

In  section  2  we  discuss  some  aspects  of  the  shear  part  of  the 
transformation  which  are  helpful  in  comprehending  the  nature  of  the  shear 
strain  that  the  transforming  particles  undergo.  In  section  3  we  propose  a 


constitutive  law  for  the  composite.  In  section  4  we  derive  a  formula  for  the 
SIF  change  due  to  the  shear  component  of  the  transformation.  Based  on  that 
formula,  we  make  first  estimates  of  the  shearing  contribution  to  SIF  change  by 
using  the  transformation  zone  derived  from  the  unperturbed  elastic  solution 
[19].  In  section  5  we  formulate  the  boundary  value  problem  for  the  stationary 
crack  and  solve  it  by  means  of  the  finite  element  method  in  section  6.  As  a 
result  the  transformation  zone  shape  and  size  are  estimated.  In  section  7  the 
estimated  zone  and  the  fracture  toughness  calculations  for  the  stationary  and 
for  the  propagating  crack  are  presented.  As  it  has  already  been  proven  in  the 
past,  the  toughness  increase  is  due  to  the  transformed  particles  left  in  the 
wake  of  the  crack  tip.  In  sections  8  and  9  the  discussion  associated  with  the 
model  results  and  the  closure  are  presented  respectively. 

2.  STRESS  INDUCED  TRANSFORMATION 

The  aim  of  this  section  is  to  introduce  basic  features  of  the  stress 
induced  transformation  so  that  the  shear  component  can  be  understood.  The 
transformation  is  martensitic  and  has  been  discussed  extensively  elsewhere 
[1,16,17,18].  It  involves  a  change  from  tetragonal  to  monoclinic  symmetry  in 
particles  in  the  composite  ceramic.  In  the  situation  of  interest  to  us,  the 
transformation  is  induced  by  critical  conditions  of  stress.  A  component  of 
the  transformation  is  a  dilatation  and  if  the  particles  were  unconstrained 
there  would  also  be  a  substantial  shear  contribution.  If  we  assume  that  the 
process  is  driven  by  the  reduction  of  free  energy  [1,7,13]  then  we  deduce  that 
the  shear  strain  component  of  the  transformation  would  align  itself  to 
maximize  the  work  done  by  the  loads  applied  to  the  particle  [7,13].  In 
addition,  the  critical  state  for  transformation  would  arise  just  when  the 
applied  loads  are  capable  of  delivering  sufficient  energy  to  the  system  to 
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compensate  for  the  increase  in  the  internal  energy  [7,13].  This  implies  a 
transformation  criterion  involving  some  combination  of  hydrostatic  and 
deviatorlc  stress.  This  issue  has  been  addressed  by  Lambropoulos  [21]. 

It  has  been  observed  that  the  situation  is  more  complicated  when  the 
transforming  particles  are  constrained  in  the  composite  matrix  [12].  The 
particles  are  capable  of  twinning  or  undergoing  some  similar  mode  of 
deformation  during  transformation.  The  twins  form  in  such  a  way  that  the 
average  shear  strain  after  transformation  can  be  quite  small  compared  with  the 
potential  unconstrained  shear  transformation  strain.  In  addition,  the  final 
strain  in  the  transformed  particle  differs  because  of  the  constraint  of  the 
surrounding  matrix.  There  is  not  yet  a  comprehensive  theory  that  accounts  for 
transformation  and  twinning.  However,  it  can  be  hypothesized  that  the 
orientation  of  the  net  shear  strain  that  results  from  transformation  and 
twinning  will  be  aligned  with  the  maximum  shear  stress  applied  by  the  matrix 
to  the  particle.  This  will  tend  to  maximize  the  external  energy  absorption 
during  transformation  and  suggests  a  critical  state  for  the  transformation 
based  on  strain  energy.  Lambropoulos  [21]  has  developed  a  constitutive  law 
for  the  constrained  transformation  along  the  above  lines  combining  the  effects 
of  hydrostatic  and  deviatorlc  stress  and  accounting  for  dilatational  and  shear 
transformation. 

In  this  paper  we  shall  use  a  simpler  law  as  a  first  step  towards  studying 
the  interactions  of  shear  transformation  with  the  crack  tip.  The  process  will 
be  considered  to  take  place  on  a  continuum  scale  and  the  description  of  the 
constitutive  law  applies  to  the  composite  ceramic.  The  criterion  for 
transformation  will  be  taken  to  be  the  achievement  of  a  critical  average 
hydrostatic  stress  in  the  composite.  This  neglects  the  contribution  made  by 
the  shear  component  to  the  work  absorbed  during  transformation.  However,  the 
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residual  shear  strains  In  constrained  transformation  particles  are  known  to  be 
small  and  the  shear  strain  In  an  unconstrained  composite  element  would  be 
correspondingly  small.  The  transformation  will  involve  a  deviatorlc  component 
as  well  as  a  dllatational  contribution.  The  deviatorlc  part  of  the 
transformation  strain  for  an  unconstrained  element  of  the  composite  ceramic 
will  be  taken  to  be  proportional  to  the  dilatational  component  of  the 
unconstrained  transformation.  All  calculations  presented  in  this  paper  will 
be  for  plane  strain  situations.  In  that  case  there  is  a  fixed  ratio  between 
the  shear  strain  and  the  volumetric  strain  in  the  transformation.  The 
orientation  of  the  shear  strain  will  be  taken  as  that  of  the  maximum  shear 
stress  when  transformation  commences  (fig.  1).  Once  the  transformation  has 
taken  place  this  orientation  will  be  locked  in  so  that  changes  of  direction  of 
the  maximum  shear  stress  will  not  cause  rotation  of  the  shear  contribution  of 
the  transformation. 

Finally,  the  transformation  will  be  supercritical  in  the  terminology  of 
Budiansky  et  al.  [20],  That  is,  at  critical  state,  the  material  transforms 
completely  without  the  existence  of  a  partially  transformed  state.  The 
details  of  the  constitutive  law  are  given  in  the  next  section. 

3.  THE  CONSTITUTIVE  LAW 

The  transformation  of  the  composite  occurs  due  to  the  martensitic 
transformation  of  the  particles.  It  takes  place  when  the  macroscopic  average 
stress  in  the  composite  is  such  that 

am  =  °mr  (1) 

m  m 

°kk 

where  *  — —  is  the  mean  stress,  cr£r  is  a  critical  value  and  o^j  is 
the  macroscopic  average  stress  tensor  at  a  point  in  the  composite. 
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The  unconstrained  transformation  strain  of  the  composite  is  which 

is  partly  dilatational  and  partly  deviatoric  in  general.  It  can  be  estimated 
to  be  the  volume  fraction  Vf  of  particles  in  the  composite  times  the 
unconstrained  transformation  strain  of  the  transforming  particles  less  the 


deviatoric  strain  nullified  by  twinning  [20],  The  supercriticality  of  the 

transformation  implies  that  there  is  no  transformation  as  long  as  am  <  a^r 

cr 

and  complete  transformation  occurs  if  a  >  a  .  The  transformation  will  be 

ram 

assumed  to  be  effectively  irreversible  in  the  conditions  prevailing  as 
observed  in  experiment  [19].  The  amount  of  the  deviatoric  strain  arising 
during  transformation  is  proportional  to  the  amount  of  dilatation.  We  will 
carry  out  calculations  for  a  variety  of  values  of  this  ratio.  The  principal 
axes  of  the  deviatoric  strain  will  be  taken  to  depend  on  the  state  of  stress 
at  transformation  as  discussed  below. 

The  stress  which  arises  [22]  in  a  constrained  element  of  the  composite  is 

°ij  =  cijki  {€kA  "  ekP 

where  is  the  final  strain  and  cijkZ  *s  tensor  of  linear  elastic 

moduli  of  the  composite  material.  For  an  isotropic  composite 

aU  '  2“  <eij  -  eIj>  +  B(ckk  -£kk)6lj  <3> 

where  p  is  the  shear  modulus,  B  is  the  bulk  modulus,  e^j  =  eij  “  y 

is  the  deviatoric  strain,  is  the  Kronecker  delta  and  eT  is  the 

dilatational  part  of  the  unconstrained  transformation  of  the  composite  (e  = 


ekk)* 


The  calculations  we  have  carried  out  are  for  plane  strain.  In  this  case 


ezz  =  0  and  we  assume  that  ez^  “  0  as  there  are  no  macroscopic  transverse 


shear  strains.  As  a  consequence 
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where  v 
written  as 


a  -  v(o 

ZZ  XX 


is  the  poisson's  ratio. 


+  a  )  -  |  p(l  +  v)  eT  (4) 

The  deviatoric  transformation  strain  can  be 


T  T  Y  .  T  Y  on 

e  =  -  e  =  -  -i—  sin2Q,  e  =  -i—  cos2u 
xx  yy  2  xy  2 

where  Q  is  the  angle  between  the  x-axis  and  a  principal  axis  of  the 
transformation  shear.  Thus,  the  transformation  is  determined  by  3  parameters 
eT,  yT  and  Q.  We  will  assume  that  Q  is  coincident  with  the  angle  to  the 
principal  axes  of  shear  stress  in  the  macroscopic  composite  state  of  stress  at 
the  instant  of  transformation  as  shown  in  fig.  1.  Calculations  are  carried 
out  for  a  variety  of  ratios  X  * 

4.  ESTIMATES  OF  THE  SHEARING  EFFECT  AT  THE  CRACK  TIP 

Before  proceeding  to  somewhat  rigorous  numerical  calculations,  we  shall 
consider  some  approximate  results  for  shearing  transformation  at  the  crack 
tip.  The  shape  of  the  zone  of  transformed  material  at  the  crack  tip  is 
determined  by  the  interaction  of  the  stresses  generated  by  the  applied  load 
and  those  generated  by  constraints  on  the  transformed  zone.  If  we  neglect  the 
latter,  we  can  estimate  the  zone  shape  as  the  locus  of  points  at  critical 
state  in  the  unperturbed  linear  elastic  solution  at  the  crack  tip.  This 
proves  to  be  quite  an  accurate  estimate  of  shape  for  the  case  of  small  scale 
dilatant  transformation  [20).  As  we  have  approximated  the  critical  state  as 
one  that  depends  only  on  hydrostatic  stress,  the  zone  shape  will  be  given  by 
the  locus  of  points  of  equal  hydrostatic  stress.  We  will  restrict  ourselves 
to  small  scale  transformation,  so  that  the  stresses  of  interest  are  given  by 
the  singular  elastic  stresses  at  the  crack  tip  due  to  the  applied  load 
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a 


ij 


KA 

—  (9) 

✓EiF  ij 


(5) 


where  KA  is  the  stress  intensity  due  to  the  applied  loads  causing  tensile 
opening  (Mode  I)  of  the  tip,  (r,9)  are  polar  coordinates  measured  from  the 
crack  tip  and  fjj  is  a  given  function  which  can  be  found  for  example  in  the 
article  by  Rice  [24],  The  shape  that  results  for  a  stationary  crack  is  shown 
In  fig.  2. 

Consider  now,  the  material  which  transforms  inside  the  zone.  If  it  were 
not  constrained  by  the  material  outside  the  zone,  a  certain  change  of  shape 
would  result.  Tractions  can  be  applied  to  the  perimeter  of  this  region  of 
material  to  return  it  exactly  to  the  shape  of  the  zone  prior  to  transfor¬ 
mation.  After  the  material  is  inserted  into  the  crack  tip  location,  the 
tractions  can  be  removed  to  give  the  final  state.  However,  the  forces 
required  to  nullify  the  constraining  tractions  T’0  will  produce  a  change  in 
the  SIF  at  the  crack  tip.  It  is  this  change  of  stress  intensity  AK  which  is 
of  interest. 

As  discussed  by  McMeeking  and  Evans  [19]  the  change  of  stress  intensity 
Is  given  by 


AK  =  /  Tc  •  h  ds  (6) 

ST 

where  _h  is  the  weight  function  [25]  whose  form  is  stated  in  McMeeking  and 

T 

Evans  [19]  and  used  with  the  assumption  that  is  homogeneous  in  A-j-  ,  the 

transforming  area  which  has  perimeter  S’j-.  We  shall  consider  now 
transformations  which  are  inhomogeneous  e^j(x).  With  the  area  A^  removed 
from  constraint  of  the  surrounding  area,  the  displacements  and  strains  [22] 

C  C 

that  result  will  be  u_  and  £  respectively  due  to  transformation  and  self- 
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sufficient  In  general  to  consider  only  the  effect  of  very  long  wakes  that 
produce  an  asymptotic  value  for  AK. 

For  X  less  than  1  the  toughness  enhancement  is  not  as  large  as  that 
given  by  equation  (18).  In  table  5  we  present  the  greatest  toughness  enhance¬ 
ments  found  for  several  X  values  and  the  amount  of  crack  advance  for  which 
they  were  obtained.  This  table  shows  that  the  toughness  enhancement  gets 
larger  when  X  gets  larger.  It  should  be  mentioned  that  the  values  for  AK 
tend  to  certain  asymptotic  values  as  Aa/w  ®  but  these  are  less  important 
than  the  maxima  presented  in  table  5.  Finally,  we  should  mention  the  face 
that  the  results  of  table  5  fit  the  curve  for  v  =  0.25  of  figure  4.  This 
means  that  the  results  taken  in  section  4  and  those  after  having  solved  the 
boundary  value  problem  of  the  stationary  crack  are  the  same. 

Applications 

Next,  we  proceed  to  see  how  the  maximum  toughness  enhancement  results 
compare  with  the  experimental  data.  We  shall  consider  v  =  0.3  even  though 
our  shearing  contribution  results  have  been  taken  by  considering  v  =  0.25. 
Figure  4  indicates  that  no  significant  difference  results.  In  the  case  of  the 
A^O^  toughened  Zirconia  whose  material  and  transformation  parameters  are 
[21]:  E  =  315  GPa,  eT  =  0.04,  w  =  10-6m,  vf  =  0.03  one  has 
i)  with  purely  dilatant  transformation 

AK  =  -  1.19  Mpa  /m 

11)  with  both  parts  of  the  transformation  and  X  =  1 

AK  =  -  1.89  MPa  /m 

ill)  Lambropoulos '  result  (both  part  of  the  transformation,  spherical 
particles  and  K^  =  5  MPa  /m) 
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based  on  equation  (16).  Such  "R  curves"  are  those  shown  in  figures  12  and  13 
for  X  »  0  and  X  =  1  respectively.  Observing  those  curves  and  the  tables  2 
through  4  we  deduce  that  the  wake  of  the  transformed  material  left  behind  the 
propagating  crack  is  a  source  of  fracture  toughness  enhancement.  First, 
figure  12  shows  that  the  toughness  enhancement  found  for  X  =  0  coincides 
with  that  found  by  McMeeking  and  Evans  [19].  ihis  coincidence  becomes  more 
pronounced  for  large  values  of  crack  advance. 

The  important  result  of  this  paper  is  the  AKg  component  that  behaves  as 
for  examaple  figure  13  shows  for  X  =  1.  As  mentioned  previously,  the  crack 
will  start  to  propagate  sooner  than  in  the  absence  of  transformation  because 
AK  is  initially  positive  due  to  AKg.  As  the  crack  grows,  AK  diminishes 
and  eventually  becomes  negative  at  about  Aa/w  =0.1.  This  means  that  the 
applied  loads  must  be  increased  to  sustain  crack  growth  and  in  a  stiff  loading 
system  the  crack  will  propagate  stably  under  rising  load.  This  will  continue 
until  Aa/w  is  about  0.7  and  the  SIF  change  is  then  given  by 


( 1— v )  AK 
X 

E  e  v^  /w 


0.35. 


(23) 


This  expression  represents  the  toughness  enhancement  and  is  larger  than  that 
given  by  equation  (18).  After  this  amount  of  crack  growth,  AK  increases  and 
so  if  the  loads  are  kept  constant  or  increased,  the  applied  K  will  exceed 
Kjq.  Thereafter  the  propagation  will  become  unstable.  Thus  the  maximum 
magnitude  of  AK  represented  by  equation  (23)  is  equivalent  to  the  asymptotic 
value  of  AK  due  to  dilatation  alone  observed  by  McMeeking  and  Evans  [19]. 

It  is  interesting  to  note  that  the  asymptotic  value  of  AK  due  to  dilatation 
plus  shear  is  not  the  relevant  quantity.  The  useful  toughness  enhancement  is 
due  to  the  minimum  in  the  curve  for  AK  at  Aa/w  =  0.7.  Thus  it  may  not  be 
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been  extended  to  the  cases  when  X  *  0  . 

In  contrast  with  AKg,  AKg  is  not  zero  and  it  does  depend  on  the  value 

of  X..  In  fact,  AKg  increases  as  X  increases.  Equation  (16)  indicates 

that  the  shearing  effect  results  in  fracture  toughness  reduction  because  AKg 

is  positive.  This  reduction  for  X  =  1  is  comparable  to  the  fracture 

toughness  enhancement  in  the  case  of  a  purely  dilatant  transformation  given  by 

equation  (18)  when  the  crack  has  advanced  Aa  *  5w  .  The  AKg  results  in 

table  1  can  almost  be  reproduced  by  the  equation  (17)  for  the  respective 

values  of  X.  This  means  that  the  transformed  material  and  the  shearing 

contribution  do  not  affect  the  features  of  the  transformation  zone  as  it  is 

found  from  the  unperturbed  elastic  crack  tip  field.  Therefore,  the  fact  that 

AKg  is  positive  is  a  consequence  of  the  nature  of  the  equation  (16).  An 

attempt  to  justify  this  may  be  made  by  regarding  the  range  of  the  positive 

contribution  to  the  integral  of  equation  (16).  Since  for  the  stationary  crack 
30 

Q  ■  —  ,  the  integral  sign  depends  on  the  sign  of  the  integrand  sin(29). 

7X  ft 

Hence,  it  is  positive  when  9  <  y  and  negative  when  9  >  y  .  Therefore,  it 
can  be  said  that  the  positive  contribution  comes  from  a  large  sector  and  can 
very  likely  override  the  negative  contribution. 

The  consequence  of  the  computed  AKg  values  would  be  that  crack 
propagation  may  take  place  sooner  than  when  there  is  no  transformation.  This 
is  because  AKg  is  positive  for  the  stationary  crack  and  K^  can  become 
greater  than  Kj£  earlier  than  otherwise.  However,  this  growth  is  likely  to 
be  stable,  a  point  which  is  elucidated  in  what  follows. 

Propagating  crack 

In  order  for  the  crack  to  advance  quasi-statically  it  is  required  that 

A  A 

the  applied  loads  K^  be  such  that  K^  =  Kj^  where  Kjq  is  given  by 
equation  (16).  Therefore,  by  knowing  AK  ,  we  can  produce  an  "R  curve"  [15] 
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zones  are  developed  inside  a  circle  centered  at  the  crack  tip  whose  area  is 
less  than  0.1 5%  of  the  mesh.  Indeed,  these  numbers  fall  within  the 
corresponding  length  range,  set  for  small  scale  yielding  condition  in  the 
elastic-plastic  fracture  analysis  by  Rice  and  Tracey  [23].  In  addition,  the 
fact  that  the  displacements  of  the  elements  near  the  perimeter  of  the  domain 
(fig.  3)  are  the  unperturbed  elastic  crack  tip  field  displacements,  ensures 
the  small  scale  condition  too.  In  figure  9,  it  can  be  seen  that  the  zone 
taken  numerically  for  X  »  0  is  almost  identical  with  that  used  by  McMeeking 
and  Evans  [19]  as  predicted  by  Budiansky  et  al.  [20].  This  means  that  the 
transformed  material  does  not  affect  the  linear  elastic  crack  tip  field  given 
by  equation  (5).  The  same  argument  applies  to  the  case  when  X  »  1  (fig. 

11).  However,  the  shearing  contribution  affects  the  zone  shape  especially  at 
small  angles  of  9.  It  is  also  worth  mentioning  that  the  zone  height  w  is 
the  same  for  all  values  of  X  as  shown  in  table  1. 

Consider  now  the  SIF  change  AKD.  It  can  be  deduced  from  table  1  that 
AKq  =0  if  we  take  into  account  the  numerical  error  involved  in  the 
calculations.  Therefore,  for  X  >  0  no  AKq  is  observed  even  though  the 
zone  boundaries  differ  from  that  with  X  =  0.  This  independence  of  AKp  from 

X  can  be  justified  by  observing  equation  (13).  The  sign  of  the  integral 

30  71 

depends  on  the  sign  of  the  cos(-2— ).  The  integral  is  positive  for  9  <  and 

negative  for  9  >  2.  .  Thus,  the  positive  contribution  to  the  integral  value 

comes  from  a  much  smaller  section  than  the  sector  of  the  negative 

contribution.  It  is  this  wider  ranging  source  of  the  negative  effect  which 

dominates  the  integral  value.  Therefore,  the  larger  range  of  r  for  small  9 

when  X  deviates  from  zero  cannot  be  cause  of  a  significant  disturbance  of 

the  final  result.  This  remark  has  already  been  made  before  [19],  but  only 


within  the  frame  of  a  purely  dilatant  transformation.  The  argument  now  has 
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Propagating  crack 

The  crack  has  been  assumed  to  propagate  quasi-statically  under  the 
conditions  mentioned  In  section  4.  The  SIF  change  calculations  AKp  and 
AKg  were  carried  out  numerically  for  a  given  crack  advance  Aa.  The  formulae 
(13)  and  (14)  were  used  again.  The  results  are  shown  In  tables  2  through  4. 

In  figures  12  and  13  the  changes  AKp  and  AKg  are  plotted  against  the  crack 
advance  Aa  for  X  *  0  and  X  =  1  respectively.  In  table  5  the  symptotlc, 
i.e.  the  maximum  SIF  change  AK  =  AKp  +  AKg  Is  shown. 

8.  DISCUSSION 

In  this  discussion  of  the  results  quoted  in  the  previous  section,  we 

shall  focus  mainly  on  how  the  shear  component  of  the  transformation  influences 

the  fracture  toughness  behavior  of  the  material.  It  should  be  borne  In  mind 

T  T 

that  the  shearing  contribution  is  characterized  by  the  parameter  X  *  y  /e  . 
Large  values  of  X  denote  large  transformation  shear  strains  for  a  given 
volume  dilatation. 

Stationary  crack 

It  has  been  mentioned  before  that  the  dilatant  part  of  the  transformation 
does  not  affect  the  material  fracture  toughness,  i.e.  AKq  *  0.  Furthermore, 
Lambropoulos  [21]  based  on  a  transformation  zone  derived  from  the  unperturbed 
linear  elastic  crack  tip  field,  concluded  that  AKg  =  0  too.  Our  analysis' 
predictions  are  in  accord  with  the  above  results  only  as  far  as  AKp  is 
concerned.  Table  1  shows  that  AKg  is  not  zero.  However,  we  have  used  a 
different  constitutive  law  from  Lambropoulos. 

Before  discussing  the  nature  of  AKD  and  AKg  we  should  emphasize  that 
the  small  scale  transformation  condition  is  satisfied.  This  is  because  the 
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Stationary  crack 


The  transformation  zone  for  X  ■  0  is  shown  in  fig.  8.  There  is  no 
shearing  effect  because  yT  -  0.  The  asterisks  indicate  the  element 
integration  stations  where  the  condition  for  transformation  was  met.  The 
transformed  zone  covers  an  area  only  0.1%  of  the  mesh.  The  nodal  displace¬ 
ments  of  the  elements  near  the  external  boundary  of  the  domain  (fig.  5)  are 
identical  within  3  significant  figures  to  the  linear  elastic  ones  when 
transformation  does  not  take  place.  The  periphery  of  the  transformation  zone 
is  shown  in  fig.  9.  The  curve  denoted  by  the  letter  D  is  the  periphery  of  the 

transformation  zone  used  by  McMeeking  and  Evans  [19].  This  last  zone  is  given 
8  2  0 

by  r  =  — —  w  cos  (•=-)  where  the  height  w  has  been  taken  to  agree  with  the 
3/3 

height  of  the  finite  element  zone. 


The  transformation  zone  for  A  =  1  is  shown  in  fig.  10.  The  boundary  of 
the  zone  is  shown  in  fig.  11  along  with  the  zone  used  by  McMeeking  and  Evans 
as  described  before.  The  zones  have  the  same  height  as  the  zone  for  A  =  0. 
However,  its  leading  front  appears  to  be  elongated  along  the  x  axis  as  y  -*•  0 
compared  to  the  other  zone.  As  a  result,  an  inflection  point  is  observed. 
The  zone  boundary  to  the  left  of  this  point  tends  to  be  concave  whereas  the 
boundary  to  the  right  tends  to  be  convex.  It  is  worth  noting  that  this 
tendency  for  elongation  increases  with  increasing  values  of  A.  This 
observation  has  been  confirmed  by  the  appearance  of  the  zone  boundaries  that 
we  generated  for  intermediate  values  of  A. 

The  SIF  change  calculations  were  carried  out  numerically  by  using  the 
formulae  (13)  and  (14)  for  the  AKD  and  AKg  changes  respectively.  These 
changes  are  shown  in  Table  1. 
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where  [ C )  Is  the  linear  elastic  stress-strain  matrix  for  plane  strain,  for  the 
composite,  {T}  is  the  two-dimensional  traction  vector,  (uN>  is  the  nodal 
displacements  and  (Et)  is  the  element  transformation  strain  vector  given  by 


2(l+v)  T 
3(l-2v)  e 


-  sin2Q 


2(l+v)  „T 

30”2v)  £ 

Y^  cos2Q 


sin2Q 


(22) 


It  is  worth  noting  that  the  stiffness  matrix  [K]  is  independent  of  the 
transformation  since  it  involves  only  the  elasticity  tensor  [C]  and  the 
interpolation  matrix  [B].  The  same  is  true  for  the  vector  {F^}  •  However, 
the  quantity  {F^}  depends  on  the  transformation  through  {E^}  . 

An  iterative  method  was  used  to  solve  equations  (21)  as  stated  in  section 
5.  All  calculations  were  carried  out  for  Poisson's  ratio  v  =  0.25.  Solutions 
were  obtained  for  \  equal  to  0,  0.25,  0.5,  0.75  and  1.  Extremely  rapid 
convergence  was  achieved  with  3  iterations  when  the  computed  displacements 
were  found  to  converge  to  within  3  significant  figures. 


7.  RESULTS 

In  this  section  we  shall  discuss  the  results  of  the  finite  element 
analysis  for  the  stationary  crack  problem.  We  shall  also  give  the  SIF  changes 
for  both  stationary  crack  and  quasi-statically  advancing  crack  deduced  through 
a  model  for  that  case.  The  shape  and  features  of  the  transformation  zone  for 
the  propagating  crack  are  determined  as  discussed  in  section  A  starting  from 
the  stationary  crack  zone.  However,  the  stationary  crack  zone  is  now  the  one 


deduced  via  the  finite  element  calculations 


(20) 


[  CkJL  Cijki  6eij  ** 
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Equation  (20)  illustrates  the  nonlinearity  of  the  problem.  That  is,  the 
second  Integral  on  the  right  hand  side  can  only  be  determined  after  the 
solution  is  found.  Thus,  the  solution  procedure  will  involve  iteration. 
Assuming  a  solution  we  determine  the  transformation  zone.  Next,  we  solve  the 
resulting  linear  problem  arising  from  equation  (20)  with  the  second  integral 
on  the  right  hand  side  evaluated  and  compare  the  solution  with  the  assumed 
one.  The  process  continues  until  convergence  is  achieved.  The  linear  elastic 
solution  can  be  used  to  initiate  the  process  with  e1  *  0  i.e.  with  no 
transformed  region.  This  method  is  addressed  in  detail  in  the  next  section. 


6.  FINITE  ELEMENT  CALCULATIONS 

The  finite  element  method  was  used  to  solve  the  boundary  value  problem  as 
stated  in  the  previous  section.  The  domain  (fig.  5)  was  discretized  into  554 
4-noded  quadrilateral  isoparametric  elements  with  4  integration  stations.  In 
figures  6  and  7  the  224  element  near  tip  mesh  and  the  320  element  far  field 
mesh  are  shown  respectively.  The  second  mesh  surrounds  the  first  one.  The 
fine  near  tip  discretization  was  used  in  order  to  ensure  that  the 
transformation  zone  is  determined  fairly  reliably.  By  using  the  standard 
element  displacement  and  strain  interpolation  matrices  [A]  and  [B] 
respectively,  we  rewrite  the  governing  equation  (20)  as  follows 


[K]  (uN)  -  (Fb>  +  (Ft) 

(21a) 

with 

[K]  -  J  [B]T  [C]  [B]  dA 

A 

(stiffness  matrix) 

(21b) 

(F  }  -  /  (A]T  (T)  ds 

D  A 

(applied  load  vector) 

(21c) 

(ft>  -  /  [b]t  {et1  ds 

(transformation  load  vector) 

( 2 1  d  ) 
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transformation  zone  appropriately*  This  is  accomplished  by  taking  into 
account  the  perturbation  of  the  linear  elastic  crack  field  Induced  by  the 
transformation.  This  subject  is  addressed  in  the  next  section. 


5.  THE  BOUNDARY  VALUE  PROBLEM  FOR  THE  STATIONARY  CRACK 

We  are  concerned  with  the  stresses  and  deformations  near  the  tip  of  a 
long  crack  in  plain  strain.  The  body  is  loaded  so  that  only  mode  I  (tensile 
opening  mode)  stresses  arise  at  the  tip.  We  are  concerned  with  “small  scale 
transformation”  where  the  transformation  zone  is  confined  to  a  region  very 
close  to  the  crack  tip.  This  can  be  achieved  by  imposing  tractions  on  a 
circular  boundary  far  from  the  tip  in  agreement  with  the  standard  singular 
linear  elastic  solution  given  by  equation  (5).  Symmetry  permits  the  analysis 
to  be  confined  to  a  semi-circular  domain  as  shown  in  fig.  5.  The  crack 
surface  is  traction  free  and  the  symmetry  line  is  free  of  shear  traction  and 
displacements  normal  to  the  line. 

The  governing  equations  in  the  plane  can  be  stated  by  the  principle  of 
virtual  work  in  the  absence  of  body  forces. 


/  C  Sty  d*  -  /  Tj  S-J  dS 

A  3 


where  A  is  the  area  and  S  is  the  perimeter  of  the  domain,  T^  the 
tractions  on  S  and  u^  the  displacements  in  A.  The  symbol  6  indicates 
an  arbitrary  virtual  variation  of  the  quantity  it  precedes.  In  addition  the 
stresses  are  related  to  the  strains  by  the  constitutive  law  given  by  equation 
(2).  As  a  result,  the  governing  equation  (19)  can  be  rewritten  as  follows 


-- '  ‘v  - 


•  V-  ^ 


.  *  .  '  .  *  .  *  fc  *  .  *  •  . 


funcdon  of  Aa  ,  the  amount  of  crack  propagation*  The  slope  versus  Aa 
decreases  continuously  and  the  curve  appears  to  approach  an  asymptote  for 
large  values  of  Aa.  At  that  stage  [19,20] 

(1-v)  AK 

- - - -  -  -  0.22  (18) 

E  e1  v^  i/w 

To  maintain  the  critical  propagation  value  of  K  ■  for  the  material  at 

the  crack  tip  the  loads  would  have  to  be  increased  with  Aa.  This  will  mean 

an  apparent  higher  value  of  K  *  K^  computed  from  the  applied  loads.  The 

material  will  appear  to  be  toughened  as  a  result. 

Consider  now  the  shearing  contribution  to  AK.  According  to  our 

hypothesis,  particles  enter  the  transformation  zone  ahead  of  the  crack  and  the 

direction  of  shearing  transformation  will  be  determined  by  the  state  of  stress 

there.  We  assume  still  that  the  stress  field  is  unperturbed  by  the  presence 
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of  a  zone  and  so  Q  equals  to  current  value  of  —  measured  from  the  crack 
tip.  As  the  crack  propagates  by,  this  value  of  Q  will  remain  unchanged  for 
that  particle.  Thus  all  particles  on  a  line  parallel  to  the  crack  will  have 
the  same  value  of  Q  ,  except  those  in  the  zone  created  before  the  crack 
propagated  (fig.  3).  The  resulting  contribution  AKg  from  such  a  wake  zone 
is  negative.  The  maximum  of  the  absolute  value  of  AKg  for  all  \  is 
achieved  when  the  crack  propagates  0.5w  whereas  the  maximum  of  the  absolute 
value  of  AKq  as  found  from  equation  (18),  is  achieved  at  Aa  *  5w.  In 
figure  A  the  ratio  |AKg  | max^  I ^olmax  *s  s^own  plotted  against  the  ratio  A. 

As  it  is  seen,  the  amount  of  the  shearing  contribution  to  the  toughening  can 
become  the  same  as  that  of  the  dilatant  contribution  when  \  approaches  1. 

In  conclusion,  the  importance  of  the  shearing  contribution  raises  the  need  for 
a  more  rigorous  estimation.  Therefore,  it  is  required  that  we  compute  the 
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where  w  is  Che  height  of  the  transformation  zone  as  shown  in  fig»  2,  X  is 

the  ratio  yT/eT  and  v  »  0.3.  The  transformation  strain  has  been  rewritten 
T  T 

38  vf  Yp  where  Yp  is  t*145  unconstrained  transformation  strain  after 
twinning  for  an  individual  particle. 

Equations  (16)  and  (17)  indicate  that  the  shearing  effect  causes  fracture 
toughness  reduction.  This  reduction  becomes  larger  as  X  increases  i.e.  the 
shear  strain  y  increases  for  a  given  eA.  As  a  consequence,  crack 
propagation  may  occur  more  readily  than  in  the  absence  of  transformation. 


Results  for  propagating  crack 

If  the  crack  propagates  quasi-statically  and  there  is  no  reverse 
transformation,  a  wake  of  transformed  material  is  left  along  the  newly  created 
crack  surfaces.  This  wake  will  be  parallel  to  the  crack  surfaces  if  we  assume 
that  the  tip  propagates  at  a  constant  value  of  K  (fig.  3).  The  dilatant 
contribution  to  AK  due  to  such  a  transformed  zone  has  been  studied  by 
McMeeking  and  Evans  [19]  and  has  been  found  to  be  a  monotonically  decreasing 


McMeeking  and  Evans  [19]  have  shown  that 


T  —  J 

AK - Li -  /  r  7  cosA  dA  (13) 

6/1 n  (1-v)  At 

where  E  is  the  Young’s  modulus.  From  the  form  of  the  weight  function  h^ 

T 

and  e  it  can  be  deduced  that 

T  -  — 

AK - ---Y . Jr2  cos(2Q  -  ^-)  sinQ  dA  (14) 

8  /lie  (1-v2)  At 
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The  fact  that  cos(— )  is  an  even  function  of  9  allows  the  integral  in 

equation  (13)  to  be  computed  in  the  region  0  <  9  <  n  (figures  2,3).  The 

same  argument  applies  to  equation  (14)  because  yT  sin9  is  an  even  function 
of  9.  That  is  because  both  y  and  sin9  are  odd  functions  of  9. 

The  SIF  change  AK  is  related  to  the  applied  field  KA  as  follows 

-  KA  +  AK  (15) 

where  Kt  is  the  S.I.F.  that  actually  describes  the  crack  tip  field  when  the 
transformation  has  taken  place.  If  AK  is  negative  equation  (15)  implies 
that  the  stress  intensity  at  the  crack  tip  is  lower  than  that  associated  with 
the  applied  loads.  In  particular,  since  the  crack  propagates  when  Kt 
becomes  equal  to  the  fracture  toughness  Kjq  of  the  p re transformed  material 
[14]  the  apparent  fracture  toughness  Kj£  is  given  by 

4  •  Kic  -  4K  <i6> 

If  AK  is  negative  the  material  appears  to  be  tougher  than  the  pretransformed 


one. 
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principle  of  virtual  work 


11C '  ‘  [  V  “9 
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Substitution  of  (8a)  into  (9)  furnishes 


ik  -J  <cL  -  eJj  -  *1?  ciju  v «  <>« 

X 

As  further  noted  by  Rice  [25]  ^i,k  a  stress  field  in  equilibrium 

in  At  because  h_  is  related  to  the  difference  between  two  displacement 
fields  each  possessing  equilibrium  stress  fields.  This  means  that  the  virtual 
work  principle  can  be  used  to  show 


/  (Eij  +  Cijkt  h*,k  “T  ■  J  “l  Cl3k*  YkK  +  “j*  dS 


R  ,  S, 


0  due  to  (8c) 


This  means  that  equation  (10)  is  written  as 


“  -  J  eO  cijki  V  “t 

*T 


For  isotropic  material,  equation  (11)  becomes 


AK  =  AKd  +  AKg 


(12a) 


where  AK^  ■  B  J  e  h^  ^  dA 
AT 


(dilatational  contribution)  (12b) 


.  »T* 

AK  ■  2p  j  e  h.  dA  (shearing  contribution) 

S  ^  a  j  1  »  J 


(12c) 


‘  ^ -  J  stl  l%V  -r .  j lVjl  1 


V  .  *■  .■ 


constraint  within  A^.  This  will  give  rise  to  an  elastic  stress  such  that 


‘ij  “  cijki  (eu  - 


aD  -  0 
ij  *1 


ni  °ij  "  °  °n  ST 


S  il 

where  n  is  the  unit  outward  normal  to  ST  and  0 . .  .  “  > -  with  x. 

—  T  ij,i  i 

denoting  position  in  a  fixed  cartesian  coordinate  system.  Tractions  are  now 

R  R 

applied  around  Aq.  producing  further  displacements  u.  anc*  8trains  £ 
that  Aj  has  its  original  shape  prior  to  transformation.  As  a  result 


+  4  -  &> 


a  *  0 

ij.i 


R  ,  S  n 
ui +  °i  ■ 0 


on  S  . 


Using  a.?  we  can  express  _T_C  as  follows 


Tc  _  _R 

T .  m  n .  g  ,  . 


i  iJ 


Thus,  equation  (6)  becomes 


“  ■  -J  "i  °ij  N 


As  noted  by  Rice  [25],  h^  can  be  treated  as  a  displacement  and  so  by  the 


iv)  experimentally  measured  value 

AK  -  -  2.4  MPa  An  (implies  X  =  1.7) 

For  another  system  whose  the  parameters  are  [21]:  E  *  470  GPa,  *  0.04 
5.10  bm,  Vf  *  0.3  one  has 

i)  with  purely  dilatant  transformation 

AK  =  -  3.96  MPa  An 

ii)  with  both  parts  of  the  transformation  and  X  *  1 

AK  =  -  6.30  MPa  An 

iii)  Larabropoulos'  result  (both  parts  of  the  transformation,  spherical 
particles  and  Kj-q  =  6  MPa  An) 

AK  *  -  4.5  MPa  An 

iv)  experimentally  measured  value 

AK  *  -  6  MPa  /m  (implies  X  -  0.9) 

In  the  case  of  the  MgO  partially  stabilized  Zirconia  whose  parameters  are 
[19]:  E  =  200  GPa,  eT  -  0.058,  w  »  6.10-bm,  vf  -  0.3  one  has: 

i)  with  purely  dilatant  transformation 

AK  =*  -0.85  MPa  An 

ii)  with  both  parts  of  the  transformation  and  X  *  1 

AK  =  -  1.35  MPa  An 

iii)  Larabropoulos'  result  (both  parts  of  the  tranf ormation ,  spherical 

*  1.97  MPa  An) 


particles  and  K 


lv)  experimentally  measured  value 

AK  *  -  2.3  MPa  /m  (implies  X  =  3) 

The  preceding  toughness  enhancement  calculations  should  be  regarded 
tentative  because  of  the  uncertainties  involved  in  the  material  parameters. 
Nevertheless,  it  is  clear  that  the  predicted  toughness  increase  may  be  in 
agreement  with  the  experimental  data  when  the  shearing  contribution  is 
considered  in  the  transformation  mechanism  analysis.  An  important  role  is 
played  by  the  parameter  X  which  determines  the  amount  of  the  shearing 
contribution.  We  have  arbitrarily  used  1  because  we  have  no  information  on 
what  values  may  be  realistic.  Alternately,  we  can  compute  a  value  of  X 
necesssary  to  bring  our  estimates  into  agreement  with  the  experimental  values. 
These  are  the  implied  values  of  X  listed.  Some  of  these  values  are  rather 
large.  However,  the  amount  of  shear  strain  during  transformation  is  typically 
quite  large  compared  to  the  dilatation.  Even  if  substantial  amounts  of  this 
are  nullified  by  twinning,  this  could  still  leave  values  of  X  of  the  order 
implied. 

9.  CLOSURE 

For  further  and  more  elaborate  treatment  of  the  shear  effect  one  may  need 
to  incorporate  size  and  particle  orientation  effects  via  a  more  realistic 
constitutive  law  for  the  composite.  This  must  include  a  consideration  of  the 
effect  of  twinning  during  the  transformation.  However,  our  simple 
calculations  have  indicated  that  shear  strain  effects  may  be  significant  and 
that  further  work  on  the  mechanics  of  the  phenomenon  may  be  profitable. 
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TABLE  CAPTIONS 


Table  1.  Stress  intensity  factor  changes  for  the  stationary  crack  calculated 
by  the  finite  element  method. 

Table  2.  Stress  intensity  factor  changes  for  the  propagating  crack  when 
X  m  y ”/e^  *  0  deduced  from  the  finite  element  calculations. 


Table  3.  Stres^  intensity  factor  changes  for  the  propagating  crack  when 
X  -  y  /e  *  0.5  deduced  from  the  finite  element  calculations. 


P* 


Table  4.  Stres^  intensity  factor  changes  for  the  propagating  crack  when 
\  =  Y  /e  =1  deduced  from  the  finite  element  calculations. 


Table  5.  Maximum  fracture  toughness  changes  achieved  after  the  crack  has 
advanced  quasi-statically. 


Table  1 


Stress  Intensity  Factor  Changes  for  the  Statlona 


(1-v)  AKjj 


(zone  height) 

T  — 

E  e1  Vj  /w 

2.14 

-0.0085 

25 

2.14 

-0.0066 

50 

2.14 

-0.0066 

,75 

2.14 

-0.0049 

2.14 

-0.0049 

Dilatational  SIF 

change 

Shearing 


SIF  change 


Table  2 


Stress  Intensity  Factor  Changes  for  the  Advancing  Crack  When  \ 


Aa 

(1-v)  AKp 

(1-v)  AKS 

(1-v)  AK 

w 

T 

E  e  /w 

E  eA  /w 

E  e  /w 

0.0 

-0.0085 

0 

-0.0085 

0.2 

-0.0780 

0 

-0.0780 

0.5 

-0.1270 

0 

-0.1270 

0.8 

-0.1550 

0 

-0.1550 

1.1 

-0.1723 

0 

-0.1723 

2.0 

-0.1961 

0 

-0.1961 

3.0 

-0.2055 

0 

-0.2055 

4.0 

-0.2097 

0 

-0.2097 

5.0 

-0.2119 

0 

-0.2119 

7.0 

-0.2142 

0 

-0.2142 

10.0 

-0.2157 

0 

-0.2157 

15.0 

-0.2167 

0 

-0.2167 

20.0 

-0.2171 

0 

-0.2171 

30.0 

-0.2174 

0 

-0.2174 

Aa 

=  Dimensionless 

crack 

:  advance 

w 

■  Dilatational 

SIF 

change 

AKS 

=  Shearing 

SIF 

change 

I 


Total 


SIF  change 


Table  3 


Stress 

Intensity  Factor  Changes  for 

the  Advancing 

Crack  When  X  *  0. 5 

Aa 

(1-v)  AKp 

(1-v)  AKg 

(1-v)  AK 

w 

E  eT  /w 

E  eT  vf  /w 

E  eT  vf  ✓ 

0.0 

-0.0066 

0.0762 

0.0695 

0.3 

-0.0961 

-0.0954 

-0.1915 

0.6 

-0.1362 

-0.1057 

-0.2419 

1.0 

-0.1656 

-0.0860 

-0.2516 

2.0 

-0.1943 

-0.0376 

-0.2319 

3.0 

-0.2037 

-0.0312 

-0.2348 

4.0 

-0.2079 

-0.0298 

-0.2377 

5.0 

-0.2101 

-0.0287 

-0.2388 

10.0 

-0.2139 

-0.0268 

-0.2407 

15.0 

-0.2149 

-0.0258 

-0.2406 

20.0 

-0.2153 

-0.0247 

-0.2400 

30.0 

-0.2156 

-0.0230 

-0.2386 

Aa 

w 

=  Dimensionless 

crack 

advance 

AKd 

=  Dilatational 

SIF 

change 

AKS 

=  Shearing 

SIF 

change 

(AK)  *  Total 


SIF  change 
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Table  4 


Stress 

Intensity  Factor  Changes 

for  the  Advancing 

Crack  When  X  ■  1 

Aa 

(1-v)  AKp 

(1-v)  AKg 

(1-v)  AK 

w 

E  e^  /w 

T 

E  e  v^  /w 

T  — 

E  e  v^  /w 

0.0 

-0.0049 

0.1527 

0.1479 

0.3 

-0.0944 

-0.1904 

-0.2848 

0.6 

-0.1344 

-0.2110 

-0.3454 

1.0 

-0.1638 

-0.1718 

-0.3356 

2.0 

-0.1925 

-0.0744 

-0.2669 

3.0 

-0.2019 

-0.0620 

-0.2638 

4.0 

-0.2061 

-0.0592 

-0.2653 

5.0 

-0.2084 

-0.0570 

-0.2653 

10.0 

-0.2121 

-0.0532 

-0.2653 

15.0 

-0.2131 

-0.0512 

-0.2643 

20.0 

-0.2135 

-0.0492 

-0.2626 

30.0 

-0.2138 

-0.0455 

-0.2594 

—  =  Dimensionless  crack  advance 
w 

AKD  *  Dilatational  SIF  change 
AKg  =  Shearing  SIF  change 
AK  *  Total  SIF  change 
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Table  5 

" Asymptotic”  Fracture  Toughness  Changes  for  the 
Quasl-Statically  Advancing  Crack 


\ 

Aa 

(1-v)  AK 

/\ 

w 

T 

E  e  /w 

0.00 

5.0 

-0.2119 

0.23 

3.2 

-0.2202 

0.50 

0.9 

-0.2527 

0.75 

0.8 

-0.2991 

1.00 

0.7 

-0.3494 

\  = 

yt/«t 

Aa 

w 

Dimensionless 

crack  advance 

AK  = 

SIF  change 

■I' ■!*.' W'pjPU.  »VAI 


i  umj  ■  pj  ■  ■:•  •■».  i  p 1 '  p  1  * 


■j  "J  ~  »■  ■  1  ».  1  vy ^ v 


t^t 


Fig.  1 

Fig.  2. 

Fig.  3. 

Fig.  4. 

Fig.  5. 

Fig.  6. 
Fig.  7. 

Fig.  8. 

Fig.  9. 
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Figure  Captions 


a.  Direction  Q  of  the  maximum  shear  stress  i 

b.  Shape  change  of  an  unconstrained  transforming  element.  The 
shearing  direction  is  along  the  maximum  shear  stress  direction. 


Transformation  zone  shape  based  on  critical  hydrostatic  stress 
transformation  and  the  unperturbed  elastic  solution  for  the 
stationary  crack.  The  maximum  height  w  of  the  zone  appears  at  4> 
=  60°.  The  shearing  direction  Q  at  transformation  is  also 
indicated  at  an  arbitrary  point  (r,0). 


Transformation  zone  shape  for  the  quasi-statically  advancing  crack. 
The  flank  of  the  zone  is  parallel  to  the  crack  surface  and  tangent 
to  the  leading  front  of  the  zone  at  9  =  60°.  The  shearing 
direction  of  the  newly  transforming  elements  is  determined  by  the 
angle  Q  at  the  leading  front  of  the  zone.  The  initially 
transformed  elements  retain  their  own  shearing  directions. 


Comparison  of  the  maximum  shearing  contribution  to  the  toughness 
enhancement  with  the  maximum  dilatant  contribution.  Those 
contributions  are  derived  from  a  zone  shape  shown  in  fig.  3  when 
the  crack  propagated.  The  zone  for  the  stationary  crack  is  shown 
in  fig.  1.  The  AKg  is  achieved  after  Aa  *  0.5w  whereas  the 
AKd  is  achieved  after  Aa  =  5w. 


Domain  and  boundary  conditions  for  the  boundary  value  problem  for 
the  stationary  crack. 


Near  tip  finite  element  mesh. 


Far  field  finite  element  mesh  which  surrounds  the  mesh  shown  in 
fig.  6. 

Transformation  zone  for  the  stationary  crack  in  the  case  of  purely 
dilatant  transformation  (y  *  0)  derived  by  solving  the  boundary 
value  problem  by  the  finite  element  method. 


Comparison  of  the  zone  boundary  of  fig.  8  with  the  zone  boundary 
derived  by  McMeeking  and  Evans  [19]  for  the  purely  dilatational 
transformation  and  marked  D.  (fig.  2). 
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Fig.  10.  Transformation  zone  for  the  stationary  crack  when  the  unconstrained 
shear  transformation  strain  is  equal  to  the  unconstrained  volume 
dilatation.  The  zone  has  been  derived  by  solving  the  boundary 
value  problem  by  the  finite  element  method. 


Fig.  11.  Comparison  of  the  zone  boundary  of  fig.  10  with  the  zone  boundary 
derived  by  McMeeking  and  Evans  [19]  for  the  purely  dilatational 
transformation  (fig.  2)  (D). 


Fig.  12.  The  R-curve  predicted  from  the  stress  intensity  factor  analysis  in 
the  case  of  purely  dilatational  transformation. 


Fig.  13.  The  R-curve  predicted  from  the  stress  intensity  factor  analysis  in 
the  case  when  the  unconstrained  shear  transformation  strain  is 
equal  to  the  unconstrained  volume  dilatation. 


SUMMARY 


A  continuum  mechanics  description,  of  the  phenomenon  of  stress  induced 
microcracking  has  been  used  to  study  the  near  tip  stress  and  strain  fields  and 
the  size  and  shape  of  a  small  scale  damaged  zone  for  a  stationary  mode  I  crack 
in  an  elastic  body.  The  material  model  is  characterized  by  a  microcracking 
criterion,  which  is  an  extension  and  simplification  of  the  generalized 
microcracking  criterion  proposed  by  Fu  and  Evans  [10]  for  the  case  of  thermal 
stress-induced  microcracking.  That  together  with  approximate  expressions 
relating  the  effective  composite  moduli  to  the  elastic  properties  of  the 
brittle  material,  via  the  microcrack  density  e  first  introduced  by  Budiansky 
and  O'Connell  [1],  yield  a  self  consistent  approach  to  the  stress  induced 
microcracking  phenomenon. 

The  raicrocrack  density  was  found  to  characterize  thr  e  regions  of 
interest.  In  the  outer  region  the  microcrack  density  is  zero  and  the  stress 
and  strain  fields  are  purely  those  for  linear  elastic  deformation.  This 
elastic  field  constrains  the  microcracking  deformation  which  in  combination 
with  material  weakening  due  to  microcracking  causes  stress  relaxation  in  a 
region  of  intermediate  raicrocracking.  Very  close  to  the  crack  tip  the 
microcrack  density  is  saturated  and  the  stress  field  becomes  again  singular 
but  with  a  lower  stress  intensity  than  would  prevail  in  the  absence  of 
microcracking.  In  the  case  where  very  rapid  microcracking  occurs  as  the 
strain  is  increased,  the  Intermediate  microcracking  zone  is  still  present 
providing  continuity  of  the  strain  field  and  a  smooth  transition  of  the  stress 
field  from  the  purely  elastic  region  to  the  region  with  saturated  microcrack 
density.  It  appears  that  the  existance  of  the  region  of  intermediate 

microcrack  density  is  essential  to  preserving  the  assumption  of  a  continuum 
composite  because  it  helps  to  avoid  any  strain  and  stress  discontinuity. 
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INTRODUCTION 


The  self  consistent  aproach  to  the  stress-induced  microcracking  of 
brittle  materials,  as  already  mentioned,  is  based  on  two  major  notions.  The 
notion  of  a  microcracking  criterion  and  that  of  effective  elastic  properties 
of  the  microcracked  material.  The  microcracking  criterion  relates  the  induced 
microcrack  density  e  to  the  existing  stress  magnitude.  The  effective  elastic 
properties  of  the  microcracked  material  are  then  functions  of  the  initial 
elastic  properties  and  the  microcrack  density.  Budansky  and  O’Connell  [1], 
introduced  the  microcrack  density  parameter  e  given  by  equation  (1) 


e 


n  P 


(1) 


where  N=nuraber  of  microcracks  per  unit  volume,  A  is  the  area  of  the  microcrack 
surface  and  P  is  the  perimeter  of  the  microcrack.  Using  energy  balance 
considerations  and  fracture  me  chan ...  s  analysis  they  were  able  to  derive 
mathematical  expressions  relating  the  effective  elastic  properties  of  the 
actual  state  to  those  of  the  unraicrocracked  elastic  state,  through  the 
microcracking  density  e.  We  will  make  use  of  the  work  of  Fu  and  Evans  [10]  to 
construct  a  microcracking  criterion  which  will  relate  Increases  in  the 

microcracking  density  to  increases  in  stress  magnitudes.  Together  with  the 

self  consistent  results  of  Budiansky  and  O'Connell  [1],  the  microcracking 

criterion  can  be  used  to  develop  a  constitutive  law  for  the  inelastic  behavior 
of  a  microcracking  continuum. 

Hur  interest  is  in  the  development  of  microcracks  near  the  tips  of  major 
cracks.  To  study  this  phenomenon  in  an  approximate  way,  we  use  the 

constitutive  law  for  the  microcracking  material  in  finite  element  calculations 
of  either  plane  strain  or  plane  stress  deformation  near  the  crack  tip.  These 
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results  are  relevant  to  the  question  of  microcrack,  toughening  in  brittle 
materials.  Before  we  proceed  to  describe  these  results  we  will  discuss 
briefly  some  previous  work  on  this  problem. 

Hoagland,  et  al.  [2],  proposed  a  simplified  method  for  estimating  the 
density  of  microcracks  in  the  vicinity  cc  a  major  crack  tip.  They  assume  a 
set  of  randomly  oriented  lines,  which  are  the  traces  of  the  potential 
microcrack  planes,  and  which  become  microcracks  when  the  stress  normal  to  the 
line  exceed  a  critical  value.  They  also  assume  that  the  dimension  of  the 
microcracks  is  very  much  less  than  the  major  crack  size  of  the  model.  In 
determining  the  stress  Intensity,  they  used  the  singular  elastic  stress  field 
near  the  crack  tip.  Later  studies  done  by  Evans  [31,  suggest  that  the  uniform 
microcrack  density  analysis  presented  by  Hoagland  et  al.  is  inaccurate  because 
there  is  an  increased  density  of  microcracks  at  the  higher  stress  levels, 
which  more  strongly  biases  the  microcrack  density  towards  the  crack  tip. 

In  1980  Hoagland  and  Embury  [4],  suggested  a  numerical  procedure  for 
modeling  the  number  and  distribution  of  microcracks  around  a  crack  tip,  as  a 
function  of  the  applied  stress  intensity.  The  procedure  accounts 
approximately  for  microcrack-microcrack  and  microcrack-crack  interactions. 
This  work  is  actually  an  extension  of  their  original  concept  of  nucleating 
microcracks  based  on  a  uniform  fracture  stress  by  adding  new  stress  fields, 
whenever  a  new  microcrack  is  nucleated.  Although  Hoagland  and  Embury  arrived 
at  some  useful  conclusions  concerning  the  stress  induced  microcracking  in  the 
vicinity  of  a  major  crack,  their  results  are  restricted  to  a  two-dimensional 
problem  with  through  thickness  microcracks  extended  along  the  third 
dimension.  In  contrast  to  that  the  theory  developed  by  Budiansky  and 
O'Connell  accounts  for  microcracks  with  finite  third  dimension  (i.e., 
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circular,  elliptical,  etc.)  and  allows  this  effect  to  be  introduced  in  an 
approximate  way. 

A  more  thorough  investigation  of  the  mechanics  of  microcrack  toughening 
was  carried  out  by  Fu  and  Evans  {10].  In  their  work  they  established  the 
basis  of  a  continuum  mechanics  description  of  the  microcrack  toughening 
process.  The  constitutive  behavior  of  the  microcrack  zone  is  based  on  the 
overall  strain  response  of  a  microcracking  medium,  subject  to  a  uniform  stress 
field.  A  generalized  microcracking  criterion  for  microcracking  at  facets 
subject  to  general  stress  is  proposed.  Stress  induced  microcracking  is 
initiated  at  a  threshold  load.  Above  the  threshold  load,  the  stress/strain 
behavior  becomes  nonlinear.  During  unloading,  linearity  resumes  and 
hysteresis  is  present,  corresponding  to  dissipation  of  strain  energy. 


THE  MICROCRACKING  CRITERION  AND  CONSTITUTIVE  LAW 

Fu  and  Evans  [10]  studied  the  phenomenon  of  stress  induced  microcracking 
of  grain  boundary  facets.  These  facets  are  generally  subject  to  residual 
stresses  caused  by  thermal  expansion  anisotropy  of  the  neighboring  grains. 
When  a  high  tensile  stress  is  applied,  favorably  oriented  facets  crack.  The 
crack  Is  confined  to  that  facet  and  it  takes  a  much  higher  stress  to  cause  it 
to  propagate  out  of  its  initial  grain  boundary  interface.  As  a  consequence, 
the  process  of  microcracking  is  at  least  initially  stable,  with  a  steadily 
increasing  stress  causing  the  nucleation  of  more  cracks  rather  than  the 
propagation  of  existing  ones. 

From  experiments  and  fracture  mechanics  analysis,  Fu  and  Evans  concluded 
that  a  particular  facet  will  crack  if  it  is  larger  than  a  critical  size  which 
depends  on  the  resolved  applied  stress  and  the  residual  stress  on  the 


interface.  In  this  notion  is  implicit  the  idea  that  facets  can  microcrack  due 


to  the  residual  stress  alone  If  the  facet  exceeds  a  critical  size.  Fu  and 

Evans  were  able  to  phrase  their  results  in  terms  of  a  parameter  £  introduced 

by  Budiansky  and  O’Connell  [1]  to  describe  microcrack  density.  When  there  are 

N  circular  raicrocracks  per  unit  volume  and  the  microcracks  are  of  uniform 

3 

size  with  radius  a  then  £  =  Na  .  Fu  and  Evans  studied  microcracking  due  to 
macroscopically  uniform  biaxial  stress  fields  (with  principal  stresses 
o^>02,  o^>0,  k=a2/0i)*  From  their  analysis  and  observations,  they  proposed 

e  =  \(o)(a1  -  oc)  for  a1  >  aQ  (2) 

where  X  is  a  parameter  that  depends  on  stress  state  and  material  properties 
but  independent  of  stress  magnitude.  For  example  X  for  k=l  is  twice  the  value 
of  X  for  k=0.  The  parameter  £  represents  the  density  of  microcracks 
introduced  by  the  applied  stress  that  is  in  excess  of  any  initial  density. 
Thus  the  microcrack  density  £  increases  linearly  with  stress  and  does  so  up  to 
a  saturation  level,  after  which  it  remains  constant.  £  is  unchanged  by  a 
reduction  of  stress. 

We  require  a  criterion  for  microcracking  which  can  be  used  in  more 
general  states  of  multiaxial  stress.  For  this  purpose  we  will  modify  equation 
(2)  to  become. 

E  “  ~  °(p  if  °R  >  aC  (3) 

where  oR  =  is  an  effective  stress,  Oq  is  the  critical  stress  for 

microcracking  initiation  and  X  is  a  material  constant  independent  of  stress 
state.  Of  course  equation  (3)  only  applies  when  the  largest  principal  stress 
is  tensile  and  increases  raonotonically .  The  criterion  (3)  is  somewhat 
different  from  that  proposed  by  Fu  and  Evans  [10].  If  we  take  the  uniaxial 
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stress  case  with  a.  *  a  as  a  datum  then  both  criteria  give  e  ■  X  (o  -  o_) 

where  X^  is  the  value  of  X  from  Fu  and  Evans  theory  for  uniaxial  tension.  For 

equal  biaxial  tension  (Cj  ■  m  0),  Fu  and  Evans'  criterion  gives 

e  *  2X  (or  -  a.)  whereas  the  one  we  have  proposed  leads  to  e  ■  X  (/2 a  -  a„)  • 
u  C  u  C' 

Thus  our  criterion  shows  a  reduced  rate  of  microcracking  in  biaxial  tension 
compared  to  the  rate  suggested  by  Fu  and  Evans  but  microcracking  starts  at  a 
lower  stress  according  to  our  criterion.  It  can  be  argued  that  the  microcrack 
density  in  biaxial  tension  would  not  be  double  that  for  uniaxial  stress,  since 
a  given  facet  cannot  be  cracked  twice.  That  is,  unless  the  families  of  facets 
cracked  by  each  applied  stress  component  are  mutually  exclusive.  However  we 
are  interested  mainly  in  a  simple  constitutive  law  with  appropriate 
characteristics  for  microcracking  to  illustrate  what  will  occur  near  major 
crack  tips.  We  will  develop  more  exact  microcracking  laws  in  future  work. 

We  summarize  our  microcracking  criterion  as  follows. 

For  monotonically  increasing  stress,  if 
0R  <  then  e  =  0  Material  remains  unmicrocracked 

°C  <  °R  <  °M  then  e  *  X(°r  ”  0q)  The  microcrack  density  is 

increasing  linearly  with  a  (4) 

°R  ^  °M  t*ien  e  =  ~  °c)  *  The  microcrack  density  is  saturated. 


Furthermore  e  cannot  decrease.  The  stress/strain  relation  given  for  a 
macroscopic  element  of  material  is  linear  for  non-mi crocracked  regions  (e»0), 
becomes  nonlinear  for  regions  with  increasing  microcrack  density  0<e<ejj  and 
again  becomes  linear  for  fully  microcracked  regions  (e  *  eM^‘  During 


unloading  linearity  resumes  and  when  material  unloads  to  zero  stress  no 
permanent  strains  are  present  as  shown  in  fig.  (1). 
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It  should  be  noted  that  the  microcrack  distribution  produced  by  a 
devlatorlc  applied  stress  coupled  with  the  residual  stress  will  be 
anisotropic.  However,  Fu  and  Evans  used  the  isotropic  theory  of  Budlansky  and 
O'Connell  [1]  to  determine  macroscopic  moduli  for  the  mlcrocracked  material, 
that  is  they  used  e  as  If  It  represented  a  random  distribution  of  microcracks 
rather  than  an  oriented  one.  We  will  follow  Fu  and  Evans  In  this  regard  to 
obtain  an  approximate  continuum  theory  for  microcracking. 

We  have  found  that  we  can  approximate  Budlansky  and  O'Connell’s  results 


as  follows. 


1-1-1  -  E  -  I 

E  v  9  f 


where  E,  v  and  E,  v  are  Young's  moduli  and  Poisson's  ratio  for  the  material 
before  and  after  stress  induced  microcracking  respectively.  Note  that  E  and 
v  are  defined  at  constant  £.  Thus  the  constitutive  equation  for  the 
microcracking  solid  becomes 


t+v  „  v  „  , 

E  ij  “  E  kk  6lj 


where  g,  Is  the  macroscopic  strain,  a  Is  the  macroscopic  stress,  Is  the 

kronecker  delta  and  rules  for  f  are  given  by  equations  (3),  (4)  and  (5). 

The  constitutive  equation  (6)  is  only  valid  for  values  of  the 
microcracking  parameter  f  within  the  Interval  i<f<“,  fig  (3).  This  together 
with  equation  (5)  defines  the  limiting  values  of  the  microcracking  density  e, 

Q 

0  <  e  <  -r-r  (See  fig  (3)).  The  upper  bound  of  e  with  the  aid  of  equation  (3) 
lo 

yield  a  relation  between  the  critical  stresses  o^,  (?M  and  the  parameter  X. 

i«e«,  ”  °C  ^  T5"  *  "JC  ^ee  ^8  (4)). 
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FORMULATION  OF  CRACK  TIP  PROBLEM 

The  high  stresses  near  the  tip  of  a  long  crack  will  cause  microcracking 
in  the  near  tip  region  in  the  material  of  interest  to  us,  as  shown  in  fig 
(5).  We  will  restrict  our  attention  to  the  situation  where  the  zone  of 
microcrack  damage  is  very  small  compared  to  the  body  containing  the  crack.  We 
will  refer  to  this  situation  as  small  scale  microcracking.  In  this  case,  the 
microcrack  zone  will  lie  within  a  region  of  uncracked  material,  some  distance 
outside  the  zone,  the  stresses  will  be  almost  the  same  as  when  the  material 
does  not  microcrack.  When  small  scale  microcracking  prevails,  these  stresses 
will  be  the  singular  linear  elastic  crack  tip  stresses. 

It  follows  that  the  plane  small  scale  raicrocracking  problem  can  be  solved 
by  considering  a  plane  region  around  the  crack  tip  to  which  are  applied 
boundary  tractions  given  by  the  linear  elastic  stress  field,  i.e. 


L  ■  "j  "u 


/ 2itr 


nj  ^ij(9) 


(7) 


where  T  are  the  boundary  tractions,  ri  is  the  outward  unit  normal  to  the 
region  boundary,  Kj  is  the  mode  I  (tensile  opening  mode)  stress  intensity 
factor  and  (r,9)  are  polar  coordinates  originating  at  the  crack  tip  as  shown 
in  fig  (6).  The  function  £  determines  the  angular  distribution  of  stress  and 
its  form  can  be  found  in  the  article  by  Rice  [23].  The  crack  surfaces  are 
traction  free.  To  ensure  small  scale  microcracking,  Kj  must  be  limited  to  a 
level  low  enough  that  the  zone  size  is  small  compared  to  the  region  for  which 
the  analysis  is  being  carried  out 

The  governing  equations  of  equilibrium  and  compatibility  are  enforced 


through  the  principle  of  virtual  work 
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u  °lj  6eij 


dA  *  /g  6u^  ds 


In  the  absence  of  body  forces,  where  A  is  the  plane  area  being  analyzed,  ST  is 
the  perimeter  where  tractions  are  prescribed,  a}  are  the  displacements  and  the 
symbol  6  indicates  a  virtual  variation  of  the  quantity  following  it  and  the 
variation  disappears  on  S  -  S-p. 

The  constitutive  law  for  the  material  has  been  described  in  the  previous 
section.  However,  we  will  find  it  useful  to  state  the  following  form. 


°ij  ci jkl  ekl  "  f+v  eij  +  (f-2v)(f+v)  ekk  6ij 


where  f  *  - j-g —  and  e  =  e(o)  as  in  equation  (3) 

1  “  T~  e 


For  plane  problems 


e  a  -  — F—  [o  fl  -  v*o  6q]  a,  p,  y  “  1»  2 
ap  E  ap  YY  ’ 


°ap  ”  f+v  leap  +  (1-iv*)  6yy 


(10a) 


(10p) 


where  v*  *  j  for  plane  strain  and  v*  =  for  plane  stress.  The  governing 
equations,  together  with  the  appropriate  boundary  conditions  corresponding  to 
a  specified  geometry,  define  a  boundary  value  problem  whose  solution  we  will 
obtain  numerically. 


THE  FINITE  ELEMENT  EQUATIONS 


The  finite  element  equations  can  be  derived  from  the  principle  of  virtual 
work  given  by  equation  (8).  Finite  element  interpolations  together  with 


.  .*•  .'-  o>  ,*.;v 
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equation  (8)  and  the  constitutive  law  (9)  and  (10.  (5),  give  rise  to  the 
nonlinear  finite  element  equations 

[k(un)Hun>  =■  {F}  (11) 

where  [k]  is  the  stiffness  matrix,  (un)  the  nodal  displacement  vector  and  {F} 
the  force  vector.  The  notation  [k(un)]  indicates  the  dependence  of  the 
stiffness  on  the  nodal  displacements  due  to  the  nonlinear  constitutive  law. 

In  the  finite  element  analysis,  a  4  noded  quadrilateral  isoparametric 

element  with  4  stations  for  the  integration  of  stiffness  was  used.  The  finite 
element  equations  (11)  were  solved  using  an  interative  method  updating  the 
stiffness  matrix  in  every  iteration.  The  microcracking  parameter  f  was  found 
to  be  the  root  of  a  sixth  order  polynomial  in  f  which  satisfies  the  conditions 
imposed  on  f  earlier  in  this  work  i.e.  l<f<°°.  The  coefficients  of  the 
polynomial  are  strain  dependent  and  had  to  be  recomputed  in  every  iteration 
(see  Appendix  II J. 

To  ensure  small  scale  microcracking,  the  outer  radius  of  the  mesh  was 

chosen  to  be  15  to  20  times  the  outer  radius  of  a  damaged  zone,  estimated  from 
the  stresses  in  the  absense  of  microcracking. 

RESULTS 

Stress-strain  fields 

Figs  (7)  to  (8)  display  the  stresses  and  strains  ahead  of  the  crack  tip 
for  various  choices  of  \  and  \  is  the  rate  of  increase  of  the  microcrack 

density  with  respect  to  the  effective  stress  and  cM  =  \(oM-oc)  is  the 
saturation  value  of  the  microcrack  density. 

In  the  unmicrocracked  region,  both  stresses  and  strains  agree  quite 

closely  with  the  solution  obtained  for  an  unmicrocracked  material.  The 

elastic  strain  and  stress  fields  constrain  the  deformation  due  to 


microcracking  in  the  damaged  zone 


The  above  constraint  together  with  the 


weakening  of  the  material  in  the  damaged  zone  due  to  microcracking,  causes 
stress  relaxation  in  the  region  of  intermediate  microcracking.  Very  close  to 
the  crack  tip  the  microcracking  density  is  saturated  and  the  stress  and  strain 
fields  become  singular  again.  We  observe  that  the  amount  of  stress  relaxation 
increases  with  X  and  it  reaches  its  maximum  level  for  X=  ®.  At  this  point  we 
notice  the  importance  of  the  presence  of  the  zone  of  intermediate 
microcracking,  for  its  presence  ensures  a  smooth  transition  for  both  the 
strain  and  stress  fields  from  the  outer  to  the  inner  asymptotic  fields.  The 
absence  of  the  above  region  would  lead  to  infinite  strain  and  stress  gradients 
along  the  boundary  of  the  fully  microcracked  material  and  the  unraicrocrack 
material  [see  Appendix  I].  We  also  notice  that  in  the  stress  relaxation 
region,  the  effective  stress  which  in  our  case  is  the  equivalent  stress 
oR  =  /o^  is  consistent  with  the  microcracking  law  eq  (3),  fig  (2)  that  we 
used.  In  the  particular  case  of  X=  ®,  oR=oc=1.0  as  we  expected.  In  this 
situation,  the  partially  microcracked  zone  is  like  a  perfectly  plastic  zone. 
The  zone  with  saturated  microcrack  density  is  dominated  by  microcracking 
deformation  with  larger  strains  and  stresses  at  lower  levels  than  would 
prevail  in  the  absence  of  microcracking. 

The  shape  and  size  of  the  transformed  zone 

A  typical  microcrack  zone  is  shown  in  fig.  (5).  Three  characteristic 
quantities  describing  the  microcrack  zones  that  we  studied  were  consistently 
computed.  x,,,  is  the  distance  along  the  x  axis  at  which  the  microcrack  zone 
boundary  crosses  the  x  axis.  hm  denotes  the  y  coordinate  of  the  farthest 
point  along  the  y  direction  to  experience  microcracking  and  @m  is  the  polar 
angle  corresponding  to  the  same  point.  The  boundary  of  the  inner  zone,  where 
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the  microcrack  density  is  saturated  is  also  characterized  by  the  corresponding 
quantities  xg,  hg  and  9g  fig.  (5). 

Fig.  (9)  to  fig  (9y)  show  how  the  above  quantities  depend  on  the  second 
critical  stress  cfy,  for  fixed  values  of  oQ  and  =  k(aM~ac)  .  The  size  of 
the  fully  microcracked  zone  is  inversely  proportional  to  the  value  of  the 
second  critical  stress  oM.  Similarly  high  saturation  values  of  the  raicrocrack 
density  yield  small  saturation  zone.  There  is  an  increase  in  the  size  of  the 
above  zone  in  both  directions  as  k  increases  or  aM  decreases.  In  no  case  does 
the  inner  zone  becomes  identical  with  the  whole  damaged  zone  because  as  we 
said  earlier  the  existences  of  a  zone  of  intermediate  microcracking  is  an 
essential  feature  of  the  microcracking  law  eq(3).  The  microcrack  zone 
boundary  depends  on  the  first  critical  stress  Oq.  For  fixed  values  of  both 
the  first  critical  stress  and  the  saturation  microcrack  density  the 

microcrack  boundary  tends  to  stretch  along  the  x  axis  while  reduction  of  1^  is 
observed,  as  k  increases.  The  angles  0m  and  9g  corresponding  to  hffl  and  hg 
drop  from  78°  for  small  values  of  k,  to  49°  as  k  becomes  infinite.  Fig  (10) 
shows  microcrack  zones  obtained  through  our  finite  element  analysis.  Fig  (11) 
and  (12)  show  the  profile  of  the  microcrack  density  ahead  of  the  crack  tip  for 
various  values  of  cfy. 

Plane  stress  results 

Similar  results  to  those  of  the  plane  strain  case  were  obtained  for  the 
plane  stress  case.  The  transverse  constraint  of  the  plane  strain  case 
produces  large  transformation  zones.  the  overall  response  of  a  microcracking 
material  was  found  to  be  similar  for  both  plane  strain  and  plane  stress  cases. 
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DISCUSSION 

Using  the  concept  of  the  microcrack  density  e  first  introduced  by 
Budiansky  and  O'Connell  [1]  together  with  a  further  simplified  microcracking 
law  proposed  by  Fu  and  Evans  [10],  we  establish  the  basis  of  a  continuum 
mechanics  description  for  the  phenomenon  of  stress-induced  microcracking  in  a 
brittle  polycrystalline  composite.  With  appropriate  integration  of  the 
constitutive  relations  we  were  able  to  obtain  information  as  to  what  happens 
in  the  process  zone  at  the  vicinity  of  the  crack  tip  of  a  major  crack  under  a 
mode  1  loading  in  a  material  susceptible  to  microcracking.  For  appropriate 
choices  of  parameters,  the  transformation  zone,  as  in  the  case  of  small  scale 
yielding  in  a  ductile  material,  is  well  contained.  The  stress  and  strain 
fields  exhibiting  the  characteristics  shown  in  fig  (7)  and  (8)  are  well 
behaved.  When  the  first  critical  stress  Oq  and  the  saturation  value 
eM  are  specified,  the  second  critical  stress  oM  or  the  variable  X  appearing  In 
the  microcracking  law  control  the  rate  at  which  the  material  microcracks  and 
affects  the  stress  and  strain  fields  as  well  as  the  shape  and  size  of  the 
transformation  zone. 

At  this  point,  we  have  all  the  necessary  tools  for  a  full  scale 
investigation  of  the  phenomenon  of  microcracking.  Such  an  investigation 
should  address  the  question  of  the  possible  effects  on  material  toughening  due 
to  microcracking.  Further  work  is  necessary  to  elucidate  the  nature  of  the 
crack  propagation  criterion.  In  addition  studies  done  by  McMeeking  and  Evans 
[6]  and  Budiansky  et  al.  [7],  suggest  that  the  toughness  of  certain  ceramics 
can  be  substantially  enhanced  through  the  controlled  use  of  martensitic 
transformation.  Faber  [9]  found  that  microcracking  contributes  to  the 
toughness  of  Zr02  ~  ceramics  and  Fu  and  Evans  [10]  found  that  stress  Induced 
microcracking  enhances  the  toughness  of  brittle  polycrystalline  aggregates. 


The  finite  element  analysis  will  substantially  contribute  to  a  better 
understanding  of  the  mechanics  of  stress  induced  microcrack  toughening  in 
brittle  material.  As  we  saw  earlier  in  our  analysis,  a  zone  with  saturated 
microcrack  density  always  exists.  We  may  be  able  to  treat  this  zone  as  being 
composed  of  homogeneous  elastic  material  and  use  a  Griffith  criterion  for 
propagation  of  the  major  crack  into  this  zone.  McMeeking,  in  unpublished 
work,  has  shown  that  the  material  obeying  the  constitutive  law  of  eq.  (3-6)  is 
nonlinear  hyperelastic.  It  follows  that  the  J-integral  is  path  independent 
throughout  the  material  in  both  unmicrocracked  and  microcracked  regions.  As  a 
consequence,  the  energy  release  rate  for  the  microcracking  material  will  have 
the  same  value  as  in  the  nonmicrocracking  material  at  the  same  applied  loads. 
Thus  if  the  critical  value  for  propagation  is  unchanged,  then  no  microcrack 
toughening  can  occur.  However,  the  critical  energy  release  rate  for  fully 
microcracked  material  may  differ  from  that  for  the  unmicrocracked  case.  A 
simple  model  of  the  crack  growing  from  microcrack  facet  to  microcrack  facet 
suggests  that  the  average  critical  energy  release  rate  would  be  less  and  so 
microcrack  embrittlement  would  occur.  However,  aspects  of  R-curve  type 
behavior  and  more  realistic  microcracktng  laws  may  hold  the  clue  to  microcrack 
toughening.  Some  of  these  issues  have  been  addressed  by  Faber  [9]  and  we  will 
consider  such  issues  in  the  future. 
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APPENDIX  I 


We  consider  an  element  consisting  of  two  dissimilar  material 


o  f  E,v 


E,  v  are  the  Young  modulus  and  Poisson's  ratio  for  an  unmicrocracked  material. 
E,  v  are  the  corresponding  moduli  for  a  microcracked  material. 

We  require  that  the  material  remains  continuous  after  deformation.  Thus  from 
equilibrium  and  compatibility  we  get  the  interface  conditions 


a  =  a 

XX  XX 


e  =  e 

yy  yy 

cr  =o  =  0 
xy  xy 


interface  conditions. 


The  constitutive  relations  for  plane  elasticity  case  are 


eafl  ■  I-  '"ae  - v*  °yy  6«P 


v*  =  v  Plane  strain 


unmicrocracked 


solid 


aaP  =  T+7  (eap  +  0-7v*T  eYY  V  v*  =  1+7  Plane  stress 


microcracked 


solid 


e«8  '  ~E~  (°ap  *  V  v*  *  T  P1“”e  3trii" 
“afs  ‘  7T7  +  ‘yy  V  7*  '  TT7  plane  8"e8s 


where  f  = 


1  ~  T~  £ 


,  where  e  is  the  microcrack  density. 


exx  strain  ahead  of  the  crack  tip 


x/( k/oc) 


Fig.  7/S  Equivalent  stress  ahead  of  the  crack  tip 


x/( k/oc) 


Fig.  7a  ayy  stress  ahead  of  the  crack  tip 


oo/XXo 


x/(. k/oc) 


ox) <  stress  ahead  of  the  crack  tip 


x/(k/oc) 


Flfl.6  The  BOUNDARY  CONDITIONS  u»od  In  solving 


th«  FINITE  ELEMENT  oquatlons 


* 

LOAD 


5  Microcracking  at  tha  vicinity  of  tha  crack  tip 


of  a  aoml-  Infinite  macrocrack  under 
mode  I  Loading 


Relations  between  the  constants  c 


<jc ,  am  and  X 


VV*.  V ' 
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APPENDIX  II 


The  raicrocracking  law  given  by  equation  (3)  earlier  in  this  paper  allow 
us  to  determine  the  microcrack  density  e  for  a  given  stress  state.  It  is  more 
convenient  to  determine  e  or  f  *  1/(1  “  -5—  e)  for  a  given  strain  state. 


E  v 

from  the  constitutive  law  {e^  +  j-jj-  e^  6^] 


we  get 


11/2._L.Ie  «  +  „  e2  1V2 

1  te.  .  e.  .  v  2  ekkJ 


o  —  [  o  o  ]  ^  ~  - -  f  e  e 

°R  1  ij  °ijJ  f+v  ieij  eij 


(f-2v) 


then  we  have  the  following  cases  in  computing  f. 


i)  aR(f)  <  ac  =>  e  =  0  => 


f  =  “  TT"  =  K 

1  -  9“e 


ii)  aR(f)  >  >  e  eM  =  \(aM  °c)  >  f  16 

1  -  9“£M 


iii)  ac  <  aR(f)  <  oM  =>  e  =  X(aR  -  aQ). 


substituting  aR  from  equation  (2)  in  eqn  (5)  we  have 


9  , .  K  »  E  f  ,  2f-v  2  }h  . 

ir <l  '  r>  ‘x  r»  leij  eij  +  v  ekt'  xv  <6> 


which  is  a  sixth  order  polynomial  in  f  of  the  form 


a6f6  +  a5f5  +  a4f4  +  a3f3  +  a2f2  +  a^  +  aQ  »  0  (7) 


whose  coefficients  are  given  below. 


stress  discontinuity 


we  allow  axx  *  oxx  and  require  eyy  -  eyy 


we  get  (f+v)(l-v*)o  =  (l-v*)o  +(v*(f+v)-v*(l+v) )o 

yy  yy  xx 


or 


and 


-  1 
a  *  7  0 

yy  f  yy 


for  Plane  stress 


yy 


,  2  2 
£  40V v 


.2  2  yy 

f  -V 


— =-  a  for  Plane  strain 

2  2  xx 


f  -v 


strain  discontinuity 


we  allow  eyy  =  eyy  and  require 


we  get 


Plane  strain: 


xx 


,  v*  (f+v)(l-2v*)  v*  1 
l!_7*  (i-^)0-2v*y-  1_7*J 


_  1-v*  (f+v )(l~2v*) 

*yy  j_-*  ( 1— v )( l-2v*)  exx 


,  f-2v  1  *  .  1-v  (f— v )(f-2v ) 

exx  =  Hl-2v  )<;  1— vT  f— v  eyy  f-v  (1— v)(  l-2v )  xx 


Plane  stress: 


“  V  ,f  “I  ^  1  f+V  V 

exx  f  j_y2  eyy  f  1+v  exx 


Notice  that  for  f=l  (unmicrocracked  material)  cry  =  cfy  and  e^  exx,  both 
stress  and  strains  are  continuous,  for  microcracked  materials  f>l  and 


a  <o  whereas  e  >  e  .  The  strains  are  larger  In  the  microcracked  region 
y  y  xx  xx 

and  the  stresses  fall  to  a  lower  value. 
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Maximum  length  of  microcraced  zone 


Maximum  length  of  saturated  zone 
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The  transformation  zone  around  the  crack  tip 


Fig.  10 $  The  transformation  zone  around  the  crack  tip 
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Fig.  lOe  The  transformation  zone  around  the  crack  tip 
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SUMMARY 

A  solution  is  presented  for  the  elastic  stress  intensity  factors  at  the 
tips  of  a  slightly  kinked,  partially  closed  crack  in  compression.  The  solu¬ 
tion  is  accurate  to  first  order  in  the  deviation  of  the  crack  surface  from  a 
straight  line  and  is  carried  out  using  perturbation  procedures  analogous  to 
those  of  Banichuk  [5],  Goldstein  and  Salganik  [6]  and  Cotterell  and  Rice  [7] 
for  the  problem  of  an  open  crack.  Comparison  with  the  exact  solution  indi¬ 
cates  that  the  asymptotic  solution  is  accurate  for  values  of  the  angle  between 
the  straight  crack  and  its  out-of-plane  kinks  up  to  about  20°. 

1.  INTRODUCTION 

Experiments  on  glass  plates  containing  pre-existing  planar  through  cracks 
oriented  at  an  angle  to  the  direction  of  the  axial  compression  have  revealed 
that  the  relative  sliding  of  the  faces  of  the  pre-existing  cracks  does  not 
result  in  co-planar  crack  growth,  but  rather  produces  at  the  tips  of  the  pre¬ 
existing  cracks  small  tension  cracks  which  deviate  at  sharp  angles  from  the 
sliding  plane  [1-4].  These  experiments  are  designed  to  be  models  for  the  pro¬ 
pagation  of  cracks  in  rocks  in  compression.  In  this  paper,  we  are  concerned 
with  the  calculation  of  stress  intensity  factors  at  the  tips  of  the  kinked 
open  extensions  of  a  closed  sliding  through  crack.  The  same  method  can  be 
extended  to  a  curved  crack  with  several  closed  sections.  The  solution  ob¬ 
tained  is  accurate  to  first  order  in  the  deviation  of  the  crack  surface  from  a 
straight  line  drawn  between  the  kink  tips  and  is  carried  out  using  perturba¬ 
tion  procedures  similar  to  those  used  in  Refs.  [5-9]  for  the  problem  of  the 
open  crack.  The  results  can  be  stated  in  terms  of  known  solutions  for  a 
single  straight  crack  or  a  co-linear  array  of  straight  cracks. 


,  AD-A154  049 


THE  MECHANICS  OF  INTERNAL  STRESS  AND  MICROCRACK 
TOUGHENING  MECHANISMS  IN  CERAMICSCU)  ILLINOIS  UNIV  AT 
URBANA  DEPT  OF.  THEORETICAL  AND  APPLIED  MECHA.  .  ,  . 


UNCLASSIrTED  */>»  C 


END 

zZC  4  v-f ,  !■:  :-u  £ 


rue v  ncunn.  .  .  y 

1/-K  0*66  ffr  />/*. 


2 


A  complete  solution  to  the  problem  of  the  sliding  kinked  crack  has  been 
given  by  Nemat-Nasser  and  Horil  (3],  who  used  a  continuous  distribution  of 
dislocations  to  model  the  crack  and  its  kinks*  In  order  to  find  the  stress 
intensity  factors,  they  solved  numerically  a  singular  integral  equation  for 
the  dislocation  distribution.  In  contrast,  we  can  avoid  the  solution  of  the 
singular  integral  equation  by  using  the  results  of  the  asymptotic  analysis  for 
the  stress  intensity  factors.  However,  the  validity  of  the  asymptotic  solu¬ 
tion  is  limited  to  small  deviations  of  the  crack  surface  from  a  straight  line. 
Comparisons  with  the  exact  solution  given  in  Ref.  [3]  indicate  that  the  first 
order  solution  for  the  mode  I  stress  intensity  factor  is  accurate  for  values 
of  the  angle  between  the  straight  crack  and  its  out-of-plane  kinks  up  to  about 
20°. 

2.  GENERAL  FORMULATION  OF  THE  PROBLEM 


2.1  Formulation  of  the  boundary  value  problem 

Consider  an  infinite  plate  of  a  homogeneous,  isotropic,  linearly  elastic, 

brittle  solid  containing  a  curved  crack  on  y  *  \(x)  ,  with  its  tips  at  the 

positions  x  =  ±a  (Fig.  1).  A  uniform  state  of  stress  a  “  ,  a  “  and  0  “ 

ax  yy  xy 

is  applied  at  infinity,  with  a  “  <  0  and  a  ®  <  0  ,  where  tension  is 

yy  xy 

regarded  as  positive.  The  corresponding  two-dimensional  boundary  value 
problem  is  given  by 

"  0 

2eij  “  “i.j  4  uj,i  in  V  *  <l> 


°ij  “  Cijk£ekJl 


°ij  "  aij 


at  infinity. 


(2) 


onn(x,  X)  ■  ong(x,  k)  *  0  on  the  open  portions  of  the  crack,  (3) 

ons(x,  X)  ■  ponn(x,  X)  on  the  sliding  portions 

u*(x,  X)  =  u”(x,  M  the  crack,  (4) 


where  Ojj  ,  e^j  and  u^  are  the  stress,  strain  and  displacement  fields  in 
the  region  V  occupied  by  the  body,  is  the  fourth  order  tensor  of  the 

elastic  moduli,  onn  and  ong  are  the  normal  and  shear  tractions  at  the 
crack  surface,  un  is  the  displacement  in  the  direction  normal  to  the  crack 
surface,  p  is  the  coefficient  of  friction,  is  QA/Sx^  and  the  super¬ 

scripts  plus  and  minus  denote  the  value  of  the  indicated  quantity  on  the  upper 
and  lower  surfaces  of  the  crack.  Note  that  the  open  and  sliding  portions  of 
the  crack  are,  in  general,  not  known  in  advance  and  their  determination  be¬ 
comes  part  of  the  solution. 

2.2  Small-parameter  expansion 

The  essence  of  the  approximation  we  use  is  that  the  solution  to  the  pro¬ 
blem  with  the  curved  crack  is  close,  in  some  sense,  to  the  solution  of  a 
similar  problem  for  a  straight  crack.  In  fact,  we  shall  use  the  solution  to 
the  following  problem,  involving  a  flat  crack,  as  the  leading  or  zeroth  order 
approximation  in  our  expansion.  Let  and  u^^  be  such  that 
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2e 
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e(0) 

ijkJTV; 


in  V' 


(5) 


at  infinity 


(6) 


(0)  00 

and  ov  '  *  a 
ana 

o(0)(x,  0)  -  o(0)(x,  0)  -  0 
yy  xy 

o(0) (x,  0)  -  pa(0)(x,  0) 

xy  yy 


u+(x,  0)  -  u^(x,  0) 


on  the  open  portions  of  the  crack,  (7) 
on  the  sliding  portions 

of  the  crack,  (8) 


where  V'  is  the  plane  with  a  straight  slit  lying  on  the  x-axis  from  -a 
to  a  .  If  the  slope  of  the  actual  crack,  \’(x)  ,  has  order  of  magnitude 
e  «  1  at  its  largest,  then  we  can  seek  a  perturbation  expansion  in  e  for 
the  solution  to  the  problem  of  the  curved  crack,  such  that 


c  -  o<0)  +  0<l> 


e  -  e 


(0)  +e(0 


(0)  .  (1) 

u  *  u  +  u 


+  0(e2 )  , 
+  0(e2 )  , 
+  0(e2 )  , 


(9) 

(10) 

(11) 


where  o<l>  .  e<1> 


.(1) 


,  ^  and  u  are  all  0(e)  compared  to  the  leading  order 


terms.  We  mention  that  X’(x)  «  0(e)  also  means  that  \(x)/a  *  0(e)  ,  be¬ 
cause  X(±a)  ■  0  .  What  remains  now  is  the  finding  of  the  equations  and  the 
boundary  conditions  governing  and  •  We  note  at  this 

stage  that  our  approach  is  identical  to  that  of  Cotterell  and  Rice  [7],  except 
that  they  addressed  the  problem  of  a  crack  open  everywhere.  Furthermore,  they 
found  their  solutions  and  expressed  their  expansions  in  terms  of 
Muskhelishvili ’s  [10]  complex  potentials.  We  prefer  to  work  in  terms  of  fun¬ 
damental  quantities,  although  it  is  entirely  possible  that  the  partially 


closed,  slightly  curved  crack  can  also  be  solved  by  a  variation  of  the  complex 


variable  treatment  of  Cotterell  and  Rice  [5J. 


We  return  now  to  the  question  of  finding  2^^  »  2^  ^  and  .  In 


order  to  find  the  equations  and  boundary  conditions  governing  o^  , 


and  u'  '  ,  we  substitute  the  expansions  (9)— (1 1)  Into  equations  (l)-(4). 


We  also  use  the  fact  that  both  X(x)  and  X'(x)  are  0(e)  to  write 


expansions  In  e  for  the  tractions  and  displacements  on  the  crack  surface  y 


X(x)  .  Using  a  tensorlal  transformation,  we  find  that  the  normal  and  shear 


tractions  on  the  actual  crack  can  be  written  as 


ann(x*  X)  •  7  [axx<x*  X)  +  °yy(x’  X)1  + 


+  —  [a  (x,  X)  -  o  (x,  X)]co820  *  a  (x,  X)sin20  , 
i  yy  xx  xy 


a  (x,  X)  -  a  (x,  X)cos20  +  i  [a  (x,  X)  -  a  (x,  X)]sin20  , 
ns  xy  ’  2  yy  xx  * 


where  9  »  X'(x)  +  0(e)  .  Then,  using  a  Maclaurln  series  expansion  in  9 


for  sin29  and  cos29  ,  we  find 


ann(x,  X)  -  ayy(x,  X)  -  2X’(x)  oxy(x,  X)  +  0(ez)  , 


a  (x,  X)  =*  a  (x,  X)  +  X'(x)[o  (x,  X)  -  oxx(x,  X)]  +  0(e2)  . 


If  we  now  write  Maclaurln  series  expansions  in  y  for  oxx  ,  a yy  and  a^y  , 


the  last  two  equations  become 


°nn(x’  X)  "  V(X’  0)  "  X(x) 


0) 


-  2X’(x)  o  (x,  0)  +  0(e)  ,  (12) 

xy 


.  ' 


do  (x,  0) 


ans(x’  X)  *  °xy(x’  0)  "  X(x) 


^ o)  - 


-  axx(x,  0)]  -r  0(e  )  , 


where  we  have  also  used  the  equilibrium  equations  5oXy/9y  *  -6oxx/9x  and 
3°yy/ay  "  -a<JXy/3x  • 

In  a  similar  way  we  can  show  that 


un(x,  X)  -  uy(x,  0)  +  X(x)eyy(x,  0)  -  \’(x)ux(x,  0)  +  0(e2) 


Using  the  expansions  (9)-(ll),  equations  (12)-(14)  can  be  written  as 

(0)  Cl)  &tJx?>(x’  0> 

ann<X-  K)  "  °yy  (X*  0)  +  °yy  (x’  0)  "  Ux)  “ ^ - 

-  2X'(x)o(0)(x,  0)  +  0(e2)  , 
xy 

,„x  /,x  do(0)(x,  0) 

°ns(x’  M  =  (X*  0)  +  '*<*•  0)  '  Mx)  - + 

+  X’CxHo^Cx,  0)  -  c<°>(x,  0)]  +  0(e2)  ,  I 

yy 


un(x,  X)  -  u<0)(x,  0)  +  u^Cx,  0)  +  X(x)e^°)(x,  0) 
-  X'(x)u£0)(x,  0)  +  0(e2)  . 


Finally,  substituting  eqns.  (9)-(U)  and  (15)-(16)  into  the  boundary  value 
problem  formulated  in  Section  3.1  (eqns.  (I )— (A) ) ,  taking  into  account  (5)-(8) 
and  separating  zero  and  first  order  terms,  we  find  that  ^  ^  and 


should  be  the  solution  to  the  following  boundary  value  problem 
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(0)  ^°xv  (0)  (I)  2 

^  -  - 2^^'(x,0)M^'<x,0H0(c'!)  , 

a  «■  g(0)(x,0)-tj(x-a)  XX.(  - -  -K)[cf(0)(x,0)-a(0)(x,0)]+o(1)(x,0)+0(e2)  . 

rcoxy  *  ox  yy  xx  xy  ’ 


Then,  the  stress  intensity  factors  can  be  calculated  as 


+  Ki^  +  *4^  +  °(e2> 


kh  ■  u«+  </2”r  «M>  ■  C  *  +  +  o(e2)  • 

r+0 


where  ,  Kj^  and  are  the  stress  intensity  factors  for 

the  zeroth  (eqns.  (5)-(8))  and  first  order  (eqns.  (18)— (22))  problems,  and 

no  0)  .  ,,  ... 

lim+[(x  -  a)3/2  ■  ^ -  +  2(x  -  a)1/2  ^(x,  0)] 

x+a  y 

-CO),  nx 

All  *  Hm+{ (x-a)3/2  ~3x~  +  (x-a)1/2[o^)(x,0)  -  a^}(x,0)]}  . 


Using  the  last  two  equations  and  a  Williams  [11]  expansion  for  the  near 
crack  tip  stress  field,  we  can  show  that 


and  K, 


3  „  e(°> 

2  "  KII  ’ 


T  -  ^0> 
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From  the  formulation  of  the  first  order  problem  (eqns.  (18)— (22)),  it  is 
clear  that  this  can  be  considered  as  the  superposition  of  the  following  two 
problems;  problem  (i)  with  a  prescribed  normal  displacment  and  zero  shear 
traction  on  the  sliding  portions  of  the  crack  and  with  the  rest  of  the  crack 
traction  free,  and  problem  (ii)  with  a  prescribed  shear  traction  and  zero 
normal  traction  everywhere  on  the  crack  face.  Thus,  knowing  the  solution  of 
the  zeroth  order  problem  (eqns.  (5)— (8))  and  having  determined  the  sliding  and 
open  portions  of  the  crack,  we  can  determine  from  the  solution  of 

problem  (i)  mentioned  above. 

As  far  as  Is  concerned,  it  is  obvious  that  only  the  prescribed 

shear  tractions  at  the  crack  surface  of  problem  (ii)  mentioned  above  that  have 
opposite  directions  on  the  upper  and  lower  surfaces  of  the  crack  have  a  non¬ 
zero  contribution  to  •  With  the  definition 

A(x)  =  j  [A(x ,  0+)  +  A(x,  0  )]  , 


CU) 

II 


is  known  (e.g. ,  [12])  to  be 


,(1) 

"II 


/na  - 


I  oLl.}(x)  ✓ 


xy 


a  +  x 
a  -  x 


dx 


(27) 


where,  according  to  (21), 


°xy>(x)  *  ^)(x)  +  X(x)  h  ^)(X)  '  ^)(x)] 

-  \'(x)[(l  +  2p2)^°)(x)  -  ^°}(x)]  . 

yy  xx 


(28) 


On  the  other  hand,  it  is  possible  that  the  stress  field  of  the  zeroth  order 

problem,  ,  has  the  characteristic  elastic  singularity  at  several 

✓  r 
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points  in  the  interval  |x|  <  a  ;  since  derivatives  of  2^^  with  respect 

-<1) 


to  x  are  involved  in  the  formula  for 


o'  '  (eqn.  (28)),  non-integrable 
xy 


singularities  will  appear  in  eqn.  (27).  To  overcome  this  difficulty,  we 

assume,  for  the  moment,  that  the  stress  components  o^^(x,  0)  are  all 

bounded  and  differentiable  with  respect  to  x  in  the  interval  |x|  <  a  ;  this 

makes  c/^(x,  0)  also  bounded  on  the  crack  face.  In  the  case  where 
xy 

o^^(x,  0)  are  singular  at  some  point  in  the  interval  jx|  <  a  ,  the  singu- 


.(0), 


larities  are  removed  by  replacing  (x,  0)  by  bounded  functions  that 

reduce  continuously  to  zero  (or  any  other  value  that  makes  a^j^x*  0) 
continuous)  over  distances  closer  than  a  small  distance  6  to  the  point  where 
the  singularities  appear.  Later  it  is  shown  that  it  is  possible  to  let  6 
tend  to  zero,  i.e.,  effectively  to  remove  the  restriction  of  bounded  and 
differentiable  o^\x,  0)  . 


(1) 


We  return  now  to  the  calculation  of  '  .  With  the  above  continuity 

assumptions  on  we  can  integrate  by  parts  eqn.  (28)  to  find 


XD 

XI 


--L  /a  -x*[(1  +  2pV°>  -^0)]  + 

..  ■  *  L  wv  L  1  wv  w  J 


/ria  -a 


yy 


yy 


XX 


+  [-*•  ' 


XX 


a  +  x 
a  -  x 


dx  - 


(29) 


-  J  (3<°>  -  )[  i  \ 1  (a)  -  $-+  (»  -  ’ft'(a)  /  i -ZJL  dx  . 

/n  -a  “  xy  7  (.  -  x)2 


It  should  be  noted  that  aX(x)  +  (a-x)X'(a)  ■  (aX(x)  +  (a-x)X'(a))  *  0 
at  x  *  a  ,  so  there  is  no  divergence  at  the  upper  limit  of  the  second  inte¬ 
gral  in  (29).  It  can  also  be  seen  that  an  integrable  singularity  can  exist  in 
—(0) 

(x)  ,  provided  it  is  not  at  x  =  ±a  ,  as  was  also  noted  in  [7].  Specifi¬ 
cally,  in  terms  of  our  earlier  discussion,  6  can  be  shrunk  to  zero  and  in 


that  limit,  the  result  of  equation  (29)  for  approaches  the  result 

obtained  by  inserting  directly  into  (29)  the  singular,  actual  a^^(x)  •  Such 
considerations,  based  essentially  on  the  fact  that  the  final  result  of  equa¬ 
tion  (29)  for  contains  o^^(x)  only  (and  not  derivatives  with  respect 

to  x),  allow  us  to  conclude  that  (29)  is  valid  for  all  integrable  a[j^x) 
(i.e.,  not  necessarily  bounded  or  continuous). 

We  mention  again  that  part  of  the  solution  to  our  problem  is  the  finding 
of  the  sliding  portions  of  the  crack.  Having  found  the  sliding  portions  of 
the  crack  and  the  solution  of  the  zeroth  order  problem,  we  can  proceed  to 
solve  the  first  order  problem  and  use  the  formulae  given  in  this  section  to 


find  the  first  order  correction  to  the  stress  intensity  factors. 


4.  THE  PROBLEM  OF  THE  KINKED  CRACK 

A  particular  case  of  the  curved  crack  is  the  kinked  crack  shown  in  Fig. 
2.  The  shape  of  the  kinked  crack  is  given  by 


mb 

b  -  a 

(x  +  a) 

for 

-a  <  x  < 

mx 

for 

|x|  <  b  , 

mb 

b  -  a 

(x  -  a) 

for 

b  <  x  <  a 

In  this  case,  a)  =  \'(a)  +  O(e^)  -  r— —  +  O(e^)  . 

n  -  a 

Following  our  previous  discussion,  we  assume  that  both  m  and  — — — 

b  -  a 

are  0(e)  ,  which  is  equivalent  to  assuming  that  \'(x)  is  0(e)  .  We 
mention  again  that  we  are  concerned  with  the  case  where  both  o  *  and  a  “ 
are  negative  and  so,  it  can  be  assumed  that  the  portion  of  the  crack  in  the 
interval  |x|  <  b  remains  closed  during  the  application  of  the  load.  Thus, 
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the  sliding  portion  of  the  crack  is  the  interval  |x|  <  b  and  the  open  por¬ 
tions  are  the  intervals  b  <  |x|  <•  a  .  Of  course,  there  is  a  possibility  that 

oo 

the  applied  stresses  and  the  orientation  of  the  kinked  crack  are  such 

that  the  whole  crack  remains  closed  and  does  not  slide.  To  check  whether  this 
happens,  we  can  solve  the  problem  assuming  that  the  crack  opens  in  the  inter¬ 
val  b  <  | x  j  <  a  ;  if  the  calculated  Kj  at  the  tips  of  the  kinks  is 
negative,  the  tip  of  the  kinks  remain  closed  and  the  assumption  that  the  crack 
opens  in  the  whole  intervals  b  <  |x|  <  a  is  in  error. 


A. 1  Solution  of  the  zeroth  order  problem 

The  zeroth  order  problem  can  be  considered  as  the  superposition  of  the 

four  problems  shown  in  Fig.  3,  where  F(x)  is  the  distribution  of  the 

0yy\x,  0)  stress  component  of  problem  no.  1.  We  note  that  for  problem  no.  4 

the  shear  stress  on  the  crack  face,  a^\x,  0)  =  pF(x)  ,  opposes  the  relative 

xy 

sliding  of  the  crack  faces.  The  quantities  of  interest  for  each  of  the  four 
problems  mentioned  above  are  given  in  the  following.  In  the  solutions  pre¬ 
sented  in  the  rest  of  this  section,  conditions  of  plane  strain  are  assumed;  in 

v 

order  to  get  the  plane  stress  solutions  we  simply  replace  v  by  ^  +  ^  . 

4.1.1  Problem  no.  1 

The  solution  to  this  problem  has  been  given  by  Erdogan  [3]  and  is  as 
follows 


a(0)(x,  0)  = 

yy 


yy 


2  E(k)  2 


,,.2  2w  2  2 ^ 

r(b  -  x  )(a  -x) 


kTFT 


-  X 


]  =  F(x) 


(30) 


a(0)(x,  0)  =  F( x )  -  a  *  , 
xx  y  v 

a(0)(x,  0)  =  0  , 
xy 
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u(0)(x,  0+)  *  u(0)(x,  0")  , 

X  X 


e(0)(x,  0+)  -  e(0)(x,  0  ) 

yy  yy 


for  x  <  b  ; 


o^^(x,  0)  =  -o 

XX  V 


for  b  <  x  <  a  , 


where  K(k)  and  E(k)  are  the  complete  elliptic  integrals  of  the  first  and 

2  2 

second  kind  respectively  and  k  =  /l  -  b  /a  . 

Also,  the  stress  intensity  factors  for  this  problem  are 


K(0)  _  °yy’/lla  r  E(k)  I 
Ki  "  k — L1  K^kT  J  * 


and  “  0  . 


Problem  no.  2  consists  of  a  plane  strain  tension.  The  solution  to  this 
'Tem  is  quite  obvious;  therefore  we  proceed  to  problem  no.  3. 

4.1.2  Problem  no.  3 

The  solution  to  this  problem  is  known  (e.g.,  [12])  to  be 


a(0)(x,  0*) 
xx 


®  2x 

xy  ^ - 2 

/a  -  x 


a(0)(x,  0*)  =  o(0)(x,  01)  =  0  , 
yy  XV 
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Fig.  1.  Infinite  plate  containing  a  curved  crack. 

Fig.  2.  Infinite  plate  containing  a  kinked  crack. 

Fig.  3.  Superposition  used  In  the  solution  of  the  zeroth  order  problem. 

Fig.  4.  Superposition  used  in  the  solution  of  the  first  order  problem. 

Fig.  5.  Infinite  plate  containing  a  kinked  crack  oriented  at  36°  to  the 
overall  compression. 

Fig.  6.  Stress  intensity  factor  at  the  tips  of  the  kinked  crack  shown  in  Fig. 
5  (p  «  0.3). 
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Finally,  using  eqns.  (46),  (47)  and  (51)  we  find  the  stress  intensity  factor 
to  be 
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♦'<*>  -a'<’>  -Trr^ryxb  (  £ «*> 


dx  ,  o 
x  -  z  X(z) 


(52) 


where  h(x)  is  the  function  determining  the  shape  of  the  wedge  of  length  2b 


2  2  2  2 

(see  Fig.  4),  X(z)  =  /(a  -  z  )(b  -  z  )  and  the  constant  C  is  determined 


from  the  equation 


i  \ 


b  J(  2  2W  2  2.  -b 

/(a  -  x  )(x  -  b  ) 


[  j  £/<.2  -  t2xb2  - 12>  jdx  - 


-  c  -  -h(b)  . 

zo  a  o 

As  discussed  in  Section  4.2,  the  shape  of  the  wedge  for  our  problem  is 
given  by 


OO  CD 


h<«>  - .  .Vi-V  l1  :.!ila.2^JL. ,  1*1  <„ . 

/a2  -  x2 


(53) 


Substituting  (53)  into  (52)  and  carrying  out  the  integrations,  we  find 


''(z)  =  Q*(z)  =  - 


ma 

2(1  -  v) 


a  -  Lto 

*y  yy 


,,  2  2.,  2  2N 

/( z  -a)(z  -b) 


T  z/z2  -  b2  -  a/a2  -  b2  ^  z2  -  z/z2  -  b2  2va2  +  b2  ^  ^ 
[v  - =5 - o -  +  - o -  "  - 5 - +  CJ 


2  2 
z  -  a 


2a 


where 


1  r  -K(k)k  +  [(1  +  2v)K(k)  -  2E(k)]k  +  2v 
K(k)  1  2k 


+  v 


a  x/x2  -  b2  -  a/a2  -  b2 

2  2 
a  -  x 


dx 


}  • 


/(a2  -  x2)(x2  -  b2) 


/■ 
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For  a  prescribed  shear  traction  ax y(x,  0)  along  the  crack  face,  it  is 
known  (e.g.,  [12])  that 


<J>’(z)  =  Q  ’  (z) 


2x  /z2  -  a2  -a 


/  a  (x,  0)  /a2  -  x2  dx  ■ 

■*  xy  *  x  -  z 


(49) 


In  our  problem 


a  (x,  0) 
xy  ’ 


HF(x) 

0 


for  | x |  <  b  , 
for  b  <  | x |  <  a  , 


(50) 


where  F(x)  is  defined  in  eqn.  (30).  Substituting  (49)  into  (50)  and 
carrying  out  the  integration  we  find 


<t> '  ( z)  =*  Q’(z) 


-T|10  {[z‘ 

2  yylL 


a2  ^-1 

3  K^kTJ 


2  2W  2  .2. 
/( z  -  a  )(z  -  b  ) 


Finally,  using  eqns.  (46)-(48)  and  the  definition 


KI  +  iKII  ”  llm+  /2ll(x  "  a)  [ayy(x’  0)  +  i0xy(x’  0)1 
x^a 

we  find  the  results  shown  in  Section  4.1.3. 


(51) 


APPENDIX  2 

The  solution  to  the  problem  of  the  opening  of  a  finite  crack  by  a  rigid 
wedge  has  been  given  by  Markuzon  [14].  In  terms  of  Muskhelishuili 's  [10] 
complex  potentials,  the  solution  is  shown  to  be 
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APPENDIX  1 


The  problem  no.  4  in  Fig.  3  is  formulated  in  terms  of  the  complex 
potentials  $  and  <l>  of  Muskhelishvili  [10].  The  stresses  and  displacements 
can  be  expressed  as 


°xx  +  °yy  “  2 +  » 


a  -a  +  2io  «•  2[z$M(z)  +  4» * ( z ) ]  , 
yy  xx  xy  * 


2G(u  +  iu  )  -  <4>(z)  -  z#'(z)  -  <Kz)  , 
x  y 


where  k  *  3  -  4v  for  plane  strain  and  k  *  (3  —  v ) /( 1  +  v)  for  plane 
stress,  the  overbar  denotes  the  complex  conjugate  and  prime  stands  for 
differentiation  with  respect  to  z  *  x  +  iy 
Introducing  the  analytic  function 


Q(z)  ■  z4»'(z)  +  <i»(z) 


eqns.  (43)-(45)  can  be  written  as 


o  +  a  -  2 [0 ' (z)  +  ♦'(«)]  , 
xx  yy 


a  -  a  +  2ia  -  2[(z  -  z)4>” (z )  -  Q(z)  -  4>’(z)]  , 
yy  xx  xy 


2G(u  +  iu  )  ~  <$(z)  -  (z  -  z)  ♦'(z)  -  Q(z)  . 
x  y 


r~ 
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be  valid,  depend  on  both  1/ c  and  6  .  But,  roughly  speaking,  the  asymptotic 
solution  is  seen  to  be  accurate  for  values  of  9  up  to  about  20°.  Unfortu¬ 
nately,  the  values  of  Kjj  for  small  values  of  9  are  not  given  in  Ref.  [3]; 
so,  comparisons  of  the  asymptotic  result  for  Kjj  with  the  exact  solution 
were  not  possible. 

6.  CLOSURE 

A  first  order  solution  has  been  obtained  for  the  stress  intensity  factors 
at  the  tips  of  the  kinked  extension  of  a  sliding  crack.  The  validity  of  the 
asymptotic  solution  is  limited  to  kinked  cracks  with  small  deviations  from 
straightness.  There  are  several  situations  where  this  deviation  is  Indeed 
small.  As  an  example,  consider  the  case  of  a  glass  plate  or  a  rock  block  con¬ 
taining  several  small  cracks  at  different  orientations.  Under  the  application 
of  a  compressive  load,  the  cracks  with  an  angle  to  the  direction  of  compres¬ 
sion,  y  ,  greater  than  Yc  *  tan  1  ~  will  remain  closed  and  only  those 
with  Y  <  Yc  can,  possibly,  slide  and  propagate.  It  is  also  known  ([l]-[4]) 
that  these  cracks  tend  to  propagate  towards  the  direction  of  compression.  So, 
if  the  coefficient  of  friction,  p  ,  is  very  high  (which  makes  Yc  small)  the 
crack  propagation  will  create  kinked  cracks  with  small  deviations  from 
straightness.  For  situations  like  these,  the  asymptotic  results  can  be  used 
to  determine  the  stress  intensity  factors  and  to  make  predictions  for  the 
direction  of  further  propagation.  In  addition,  fatigue  due  to  non¬ 
proportional  loads  can  cause  the  development  of  cracks  that  are  not  straight 
and  are  partially  closed,  although  open  at  the  tip.  For  cases  where  the 
deviation  from  the  straight  line  is  small,  the  methods  devised  here  can  be 
used,  although  a  criterion  for  determining  where  the  closed  portions  lie  would 


have  to  be  developed 


Kn  ”Ki?)  +  7“  Ki°)  +  *11*  +  0(e2>  * 


where  and  are  given  in  equations  (38),  (39), 

(41)  and  (42)  respectively.  Equations  (41)  and  (42)  show  that  ^  and 
depend  on  Poisson's  ratio  v  .  This  means  that  they  have  different 
values  under  plane  strain  or  plane  stress  conditions.  On  the  other  hand,  it 
is  known  that,  in  the  absence  of  body  forces,  the  stress  intensity  factors  for 
traction  boundary  value  problems  are  the  same  under  plane  strain  or  plane 
stress  conditions  and  independent  of  the  elastic  constants.  The  reason  that 
v  enters  the  expressions  of  the  asymptotic  solution  for  the  stress  intensity 
factors  is  that  the  displacement  field  of  the  zeroth  order  problem,  which 
depends  on  v  ,  was  used  to  formulate  the  first  order  problem.  As  a  result, 
the  in-plane  components  of  the  first  order  correction  to  the  stress  field  and 
the  corresponding  corrections  to  the  stress  intensity  factors  depend  on  v  . 
However,  numerical  calculations  of  and  show  that  their  depen¬ 

dence  on  v  is  very  weak  and  that  their  values  for  plane  strain  and  plane 
stress  conditions  are  practically  indistinguishable,  which  validates  the  first 
order  correction. 

Next,  we  apply  our  results  to  the  problem  of  a  infinite  plate  containing 
a  kinked  crack  oriented  at  36°  to  the  overall  compression  (Fig.  5).  The  exact 
solution,  given  in  Ref.  [3],  and  the  asymptotic  results  for  Kj  are  plotted 
in  Fig.  6  versus  the  angle  between  the  straight  crack  and  its  out-of-plane 
kinks,  9  ,  for  several  values  of  the  ratio  of  the  length  of  the  kink,  A  ,  to 
the  length  of  the  straight  crack,  c  .  In  general,  the  region  of  accuracy  of 
the  asymptotic  solution  depends  on  both  i/c  and  9  ,  because  the  values  of 
ra  and  u>  in  our  analysis,  which  must  be  small  for  the  asymptotic  solution  to 


K(k)k  +  [(1  +  2v )K(k)  -  2E(k)]k  +  2v 


We  proceed  now  to  the  calculation  of  the  mode  II  stress  intensity  factor 

for  the  first  order  problem,  ,  which  can  be  determined  either  by  solving 

problem  no.  2  in  Fig.  4  or,  equivalently,  using  eqn.  (29).  We  note  that 

cr^^(x)  ,  ^(x)  and  cjl^(x)  all  have  the  characteristic  —  singularity 

xx  yy  xy  . 

n  y  r 

at  x  =  ±b  ;  but  as  discussed  in  Section  3,  eqn.  (29)  for  can  still  be 

used,  provided  that  the  singularities  are  integrable,  which  is  indeed  the 
case.  So,  using  the  solution  of  the  zeroth  order  problem  derived  in  the  pre¬ 
vious  section  and  applying  eqn.  (29),  after  some  lengthy,  but  straightforward, 
integrations  we  find 

K  =  pK^  +  (o  -  o  )  -t-— ■ —  (l  -  sin  1  —  )  +  2p2mo  /na  + 
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>re  Kj  is  given  in  eqn.  (40). 


5.  DISCUSSION  AND  COMPARISON  WITH  THE  EXACT  SOLUTION 


The  obtained  asymptotic  solution  for  the  stress  intensity  factors  at  the 
tips  of  a  kinked  crack  is  of  the  form 


K,  -  K(t0)  -  -I  a)  +  0(e2)  , 
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for  b  <  | x |  <  a  . 

Also,  the  stress  intensity  factors  for  the  zeroth  order  problem  are  given 
by 
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(38) 


pa  ")  /Ha  +  pK<0)  . 

yy  i 


(39) 


4.2  Stress  intensity  factors  for  the  first  order  problem 

As  discussed  in  Section  3,  the  first  order  problem  can  be  considered  to 
be  the  superposition  of  the  two  problems  shown  in  Fig.  4. 

The  mode  I  stress  intensity  factor  Is  determined  by  solving 

problem  no.  1  in  Fig.  4,  which  is  actually  the  problem  of  the  opening  of  a 
finite  crack  by  a  rigid  wedge.  The  general  solution  to  this  problem  has  been 
given  by  Markuzon  [14].  Taking  into  account  eqn.  (22)  and  the  solution  of  the 
zeroth  order  problem  derived  in  the  previous  section,  we  find  that,  for  our 
particular  case,  the  shape  of  the  wedge,  h(x)  ,  is  given  by 


2h(x)  =  Uy1)(x,0+)-u£1)(x,0') 


2m 
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Using  the  above  formula  for  the  shape  of  the  wedge  and  Markuzon's  [14] 
solution,  we  find  the  mode  I  stress  intensity  factor  to  be  (see  Appendix  2) 
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for  b  <  x  <  a  . 


Also, 
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4.1.4  Superposition 


Superimposing  the  solutions  of  the  four  problems  shown  in  Fig.  3,  we  find 


the  quantities  of  interest  of  the  solution  to  the  zeroth  order  problem  to  be 


_ ( CO  «*>  00 

O'  vx )  -  F(x)  -  O  +0  , 

xx  yy  xx  * 


o(0)  =  F(x)  , 

yy 


^(x)  =  pF(x)  , 


i0)(x.  0+)  -  .«»U,  0-)  .  2v 


,  2  2 
r  a  -  x 


0  -  pO 


u^(x,  0+)  -  u^°\x,  0  )  =  2(1  -  v) 
x  x  0> 


2  2 


/a  -  x 


for  x  <  b  ;  and 
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and  u^^(x,  0*)  *  ±(1  -  v)  -221  /a^  -  x^  for  x  <  a  , 
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where  G  is  the  shear  modulus  of  the  material. 
In  addition. 


K^0)  =  0  , 


and  Kj?)  =  o  /Ha 
II  xy 


4.1.3  Problem  no.  4 


The  solution  to  this  problem  is  derived  in  Appendix  1  and  is  as  follows 
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